Loop Braiding Statistics and Interacting Fermionic Symmetry-Protected Topological Phases in Three Dimensions

We study Abelian braiding statistics of loop excitations in three-dimensional (3D) gauge theories with fermionic particles and the closely related problem of classifying 3D fermionic symmetry-protected topological (FSPT) phases with unitary symmetries. It is known that the two problems are related by turning FSPT phases into gauge theories through gauging the global symmetry of the former. We show that there exist certain types of Abelian loop braiding statistics that are allowed only in the the presence of fermionic particles, which correspond to 3D"intrinsic"FSPT phases, i.e., those that do not stem from bosonic SPT phases. While such intrinsic FSPT phases are ubiquitous in 2D systems and in 3D systems with anti-unitary symmetries, their existence in 3D systems with unitary symmetries was not confirmed previously due to the fact that strong interaction is necessary to realize them. We show that the simplest unitary symmetry to support 3D intrinsic FSPT phases is $\mathbb{Z}_2\times\mathbb{Z}_4$. To establish the results, we first derive a complete set of physical constraints on Abelian loop braiding statistics. Solving the constraints, we obtain all possible Abelian loop braiding statistics in 3D gauge theories, including those that correspond to intrinsic FSPT phases. Then, we construct exactly soluble state-sum models to realize the loop braiding statistics. These state-sum models generalize the well-known Crane-Yetter and Dijkgraaf-Witten models.

Topological phases in three spatial dimensions (3D) can support particle and loop excitations 1 . While we learn in undergraduate quantum mechanics that there are only bosonic and fermionic exchange statistics for particles in 3D, the rich statistical properties of loop excitations only begin to be uncovered recently 2,3 , in conjunction with the study of bosonic symmetry-protected topological (BSPT) phases in 3D 4-6 . In particular, it was found that to characterize BSPT phases protected by finite unitary symmetries, it is necessary to consider the so-called "three-loop" braiding statistics in gauge theories, a new type of braiding statistics in 3D 2,3,7-9 , where the gauge theories are obtained by gauging the global symmetries of BSPT phases.
So far, all studies of loop braiding statistics have focused on gauge theories where all particles are bosons. Loop braiding statistics in the presence of fermionic particles are much less explored. Perhaps the most important question to address is: Does the presence of fermions allow new types of loop braiding statistics that are not possible otherwise?
This question is closely related to the problem of classifying interacting fermionic symmetry-protected topological phases (FSPT) in 3D, in which the braiding statistics of vortex loops serves as a topological invariant for the bulk phase. To put it into context, we briefly review the classification of FSPT phases with unitary symmetries in 3D. For non-interacting fermionic systems, it is well known that there are no nontrivial FSPT phases protected by on-site unitary symmetries  . On the other hand, we can create interacting FSPT phases by effectively turning fermions into bosons with the help of strong interactions (i.e. fermions forming spins or molecular bound states), and letting the bosons form SPT states. An interesting question then arises: are there "intrinsic" FSPT phases in 3D protected by unitary symmetries? Here, by "intrinsic", we mean those FSPT phases that do not stem from BSPT phases. If there were any, we know that they must be strongly interacting states (since we know there are no nontrivial non-interacting SPT phases). Then, the loop braiding statistics obtained by gauging symmetries of these "intrinsic" FSPT phases are the ones that are only allowed in the presence of fermionic particles.
Due to the non-existence of non-interacting phases, one has to confront the complexity of interacting systems from the very beginning to realize FSPT phases. A fruitful approach is to construct and study exactly-solvable lattice models. A systematic construction of fermionic SPT phases has been proposed in Ref. 14, although no explicit examples in 3D for unitary symmetries (except the bosonic ones) were given. In addition, it is not clear what kind of physical properties characterize the constructed FSPT states in Ref. 14. We should mention that other classification schemes have been proposed, such as spin cobordism group in Ref. 15

and the invertible TQFTs in
Refs. 16 and 17. In this work, we study Abelian loop braiding statistics in 3D gauge theories with fermionic particles. As we will see, the presence of fermionic particles indeed enables new types of three-loop braiding statistics, forbidden when the particles are bosonic. Because of the well-established correspondence between SPT phases and gauge theories 2,7,[18][19][20] , our results also imply existence of intrinsic FSPT phases. We derive these results through a combination of physical arguments and exactly-solvable models. Moreover, we derive a complete classification of Abelian three-loop braiding statistics in Abelian gauge theories (i.e., the gauge group is Abelian) in the presence of fermionic particles.
We show that the simplest symmetry group which allows for intrinsic interacting FSPT phases is Z f 2 ×Z 2 ×Z 4 . Here Z f 2 is the fermion parity conservation. It turns out that this example captures the essence of all 3D "Abelian" FSPT phases discussed in this work. Hence, we give an intuitive picture for one of the Z f 2 × Z 2 × Z 4 intrinsic FSPT phases -namely, the "root" phase -in terms of decorated domain walls. In this description, a symmetric state can be obtained by proliferating domain walls of the global symmetry. If domain walls themselves are "decorated" by lower-dimensional SPT phases, then the wavefunction of proliferated domain walls may represent a nontrivial SPT state 21 . In the Z f 2 × Z 2 × Z 4 FSPT phase, a Z 4 domain wall is decorated by a 2D FSPT phase protected by the Z f 2 × Z 2 symmetry. Using the terminology of Ref. 22, the one that we use for decoration is the "root" Abelian FSPT phase in 2D, which has a Z 4 classification (note that the full classification of 2D Z f 2 × Z 2 FSPT phases is Z 8 ; however, the Z 8 root phase is "non-Abelian"). In fact, because this Abelian root phase has a Z 4 classification, it can only exist on a Z 4 domain walls 14,20,22,23 . This 2D phase can be easily realized with non-interacting fermions, and a simple example is the following: the system consists of two layers, a Chern insulator with Chern number C = 1 (or equivalently, two copies of p x + ip y superconductors), and its time-reversal image with C = −1. This non-interacting FSPT has helical Dirac fermions on the edge, and since the two chiral modes carry opposite Z 2 charges, they can not backscatter. Interactions on the edge can cause spontaneous breaking of the Z 2 symmetry, but the edge can not be both symmetric and non-degenerate 20,24 .
To expose the physics of the 3D FSPT phase constructed above, we imagine inserting a Z 4 symmetry defect loop into the FSPT phase (Fig. 1). Since the defect loop can be thought of as the effective boundary of a domain wall, it must carry similar gapless modes as the edge of the 2D FSPT phase decorated on the domain wall, assuming no spontaneous symmetry breaking along the loop. As we already mentioned, one possible effective low-energy theory for the gapless modes is a 1D helical Dirac fermion. This 1D hilical Dirac fermion on Z 4 symmetry defects is an important property of the FSPT phase.
To explicitly show that the above 3D FSPT phase is intrinsic, we need to gauge the Z f 2 × Z 2 × Z 4 symmetry. Then, symmetry defect loops turn into dynamical vortex loops in the gauged system. There are two "three-loop" braiding process (see Fig. 2) which can reveal the intrinsic nature of the FSPT phase: first, consider braiding a Z f 2 (i.e. fermion parity) vortex loop around a Z 2 vortex loop, while both are linked to a unit Z 4 vortex loop. We find that this three-loop braiding statistics is either semionic or anti-semionic (±π/2). On the other hand, using a result that will be established in Sec. III D, if this FSPT phase stems from a BSPT phase, this threeloop braiding phase can only be 0 or π. The essence in this difference is that fermion parity vortex loops play a nontrivial role in the three-loop braiding statistics in the FSPT phase constructed from decorated domain walls. Hence, this FSPT must be intrinsically fermionic and the corresponding loop braiding statistics can exist only in the presence of fermionic particles. The other process is to exchange two identical Z 2 vortex loops linked to the unit Z 4 vortex. The resulting Berry phase turns out to be π/4 (up to a π ambiguity), which as we will see is not allowed in systems with bosonic charges.
Due to the length of the paper, we briefly outline the strategy underlying our approach: in Sec. III we define a set of topological invariants to unambiguously characterize Abelian three-loop braiding statistics in the presence of fermionic particles, and derive physical constraints satisfied by the topological invariants. We then solve the constraints to obtain possible solutions corresponding to gauged intrinsic FSPT phases. Next, to show that the constraints are complete and the solutions we found are indeed physical, we consider an exactly-solvable lattice model of topological twisted gauge theories which realize all solutions we have found in Sec. IV and V.

II. PRELIMINARIES
A. Symmetry Any fermionic system has a fundamental unbreakable symmetry, namely the conservation of the total fermion parity: P f = (−1) N f , where N f is the number of fermions. The two operators {1, P f } form a symmetry group, which we denote as Z f 2 . In addition, the system may be symmetric under other global symmetry transformations. Together with the fermion parity, all symmetry transformations form a symmetry group G. It is generally required that P f commutes with all elements in G. Accordingly, Z f 2 is a normal subgroup of G. The quotient group G = G/Z f 2 in a sense contains all the "physical" symmetries (i.e. those can be broken by physical perturbations). Mathematically, G is a central extension of G by Z f 2 . Given G, such an extension is not unique. Possible extensions are mathematically classified by the second group cohomology H 2 [G, Z f 2 ]. For example, it is well known that fermionic systems with time-reversal symmetry Z T 2 = {1, T } come in two types: one with T 2 = 1, and the other with T 2 = P f , corresponding to the two elements in In this work, we consider fermionic systems with an Abelian unitary symmetry group G. Without loss of generality, we can represent G as follows: 25 where N 0 is an even number. We use integer vectors a = (a 0 , a 1 , . . . , a K ) to denote the group elements, with a i = 0, 1, . . . , (N i − 1), and use additive convention for group multiplication. Generators of the group are denoted as e i = (0, . . . , 1, . . . , 0) where the ith entry equals 1 and other entries equals 0.
The group Z f N0 is singled out because the fermion parity element is given by Equivalently, it means that the unit charge under Z f N0 is a fermion, while the unit charge under other Z Ni subgroups are all bosons. If we consider the more "physical" global symmetry group G = Z N0/2 × K i=1 Z Ni , we find that all fermions carry "half charges" under Z N0/2 , forming the so-called "projective representations" of Z N0/2 . (If N 0 /2 is odd, this actually does not give a true projective representation, due to the familiar isomorphism that

B. FSPT phases and gauging symmetry
As discussed in the introduction, fermionic systems with a symmetry G may form various SPT phases, i.e., gapped symmetric short-range entangled states. In this work we study FSPT phases protected by Abelian unitary symmetries in Eq. (1).
One way to characterize FSPT phases is to "gauge" the global symmetry G in lattice Hamiltonians. That is, we minimally couple the system to a lattice gauge field of a (discrete) gauge group G. There is a well-defined procedure to do so; see e.g. Refs. 7 and 18. The resulting model, a gauge theory coupled to fermionic matter, is actually topologically ordered, in the sense that it is gapped and it hosts topologically nontrivial excitations. In a gauge theory, nontrivial excitations carry either gauge charge or gauge flux. It has been proposed and verified in various systems that braiding statistics of charge and flux excitations in the gauged models are able to distinguish the original SPT phases. In this work, we study braiding statistics in 3D gauged FSPT systems, extending previous works on 2D/3D BSPT phases 2,7 and 2D FSPT phases 22,25 .
It is sometimes useful to only gauge the fermion parity symmetry subgroup Z f N0 . In this way, the fermionic system is turned into a bosonic one, in the sense that there are no local fermionic excitations. Because of the direct product structure in Eq. (1), the other global symmetries remain unaffected by the gauging procedure, so we end up with a Z N0 gauge theory enriched by a symmetry group G/Z N0 = K i=1 Z Ni 9,26-30 .

C. Basics of 3D braiding statistics
We now discuss the basics of braiding statistics between excitations in a gauged 3D FSPT system, i.e., a G gauge theory coupled to fermionic matter.
There are two kinds of excitations in the system: particle excitations that carry gauge charges, and vortex loop excitations that carry gauge fluxes. For an Abelian group G, we use integer vectors q = (q 0 , q 1 , . . . , q K ), with q i = 0, . . . , (N i − 1), to denote the charge excitations. The self statistics associated with exchanging two identical charges is given by That is, it is a fermion/boson if q 0 is odd/even. Vortex excitations are string-like and must form closed loops inside the bulk of the system. They carry gauge flux. We use vectors to label gauge fluxes, where a i = 0, 1, . . . , (N i − 1) is an integer. There is a well-known correspondence between gauge fluxes and group elements: one may regard the vector a = (a 0 , . . . , a K ) that labels φ as an group element of G. Accordingly, the fermion parity group element corresponds to the fermion parity flux (π, 0, . . . , 0). Unlike charge excitations, vortex excitations cannot be uniquely labeled by their gauge fluxes. Two distinct vortices α and α ′ may carry the same gauge flux, i.e., φ α = φ α ′ . It can be shown that two vortices carrying the same gauge flux can be transformed to one another by attaching gauge charges. The mutual braiding statistics between a charge excitation q and a vortex loop α follows the Aharonov-Bohm law where "·" is the vector inner product. In particular, the Aharonov-Bohm phase between q and a fermion parity vortex is given by πq 0 . The most interesting part of 3D braiding statistics is between vortex loops. It was shown 2,3 that the fundamental loop braiding process involves three loops (Fig. 2): a loop α is braided around a loop β while both are linked to a third "base" loop γ. This three-loop braiding process has been used to characterize various 3D topological phases. 7 Following the notation in Ref. 2, we denote the three-loop braiding phase by θ αβ,c , where c is the integer vector that labels the gauge fluxes carried by γ: c ≡ φ γ = ( 2π N1 c 1 , . . . , 2π NK c K ). We use θ αβ,c instead of θ αβ,γ , because the three-loop braiding statistical phase depends only on the flux φ γ parametrized by c and is totally insensitive to the amount of charges attached to γ. We will also consider an exchange or half-braiding process: two identical loops α, both linked to the base loop γ, exchange their positions. We denote this three-loop exchange statistics by θ α,c .
In most part of the paper, we will only consider Abelian braiding statistics. Non-Abelian loop braiding statistics can also appear in gauge theories with Abelian gauge group. We will briefly touch upon non-Abelian loop braiding at the end of the paper.

III. PHYSICAL CONSTRAINTS ON ABELIAN LOOP BRAIDING STATISTICS
In this section, we study general properties of loop braiding statistics in gauged 3D FSPT systems. For sim-plicity, we only consider Abelian loop statistics, i.e., every Berry phase associated with braiding excitations is Abelian. We will discuss physical constraints on Abelian loop braiding statistics, and will discuss which types can result from "intrinsic" FSPT phases.

A. Topological invariants
To begin, we define a set of topological invariants {Θ ij,k , Θ i,k } for loop braiding statistics. These topological invariants are a subset of the full braiding statistics data, and hence are simpler to deal with compared to the latter. Nevertheless, they are equivalent to the full set of loop braiding statistics, since the latter can be reconstructed out of the former. Similar topological invariants have been introduced in 2D/3D gauged BSPT phases 2,7 , as well as 2D gauged FSPT phases. 25 Let α, β and γ be vortex loops, carrying unit flux 2π Ni e i , 2π Nj e j and 2π N k e k , respectively. Here, e i is an integer vector (0, . . . , 1, . . . , 0) where the ith entry is 1 and all other entries are 0, with i = 0, 1, . . . , K. We define the topological invariant Θ ij,k as follows: where θ αβ,e k is the mutual braiding statistics between α and β while both are linked to the base loop γ (Fig. 2). Here, we use N ij to denote the least common multiple of N i and N j . Similarly, we define a topological invariant Θ i,k for the self statistics associated with exchanging two identical α's, both of which are linked to the base loop γ. It is defined as follows and The above topological invariants {Θ ij,k , Θ i,k } are defined in a way such that (i) they only depend on the flux of α, β and γ and (ii) the full set of braiding statistics can be reconstructed out of {Θ ij,k , Θ i,k }. One can check the property (i) by replacing α, β, γ with α ′ , β ′ , γ ′ respectively. It is easy to show that the topological invariants do not change if φ α ′ = φ α , φ β ′ = φ β and φ γ ′ = φ γ . The proof of property (ii) is more involved, so we give the proof in Appendix A.

B. Physical constraints
The topological invariants {Θ ij,k , Θ i,k } cannot take arbitrary values. We argue that the topological invariants have to satisfy the following constraints: where N i...k denotes the greatest common divisor of N i , . . . , N k , and N i...k denotes the least common multiple of N i , . . . , N k . These constraints can be proved by checking various consistency conditions on the three-loop braiding statistics. The proofs are given separately in Appendix B. These constraints are necessary conditions that physical Abelian three-loop braiding must satisfy (except Eq. (8h) which remains a conjecture at this stage). On the other hand, we do not know at this point whether these constraints are also sufficient, in the sense that every solution to these constraints can be realized in physical systems. To verify the completeness of the constraints, we will present a family of exactly solvable lattice models in Sec. IV and show that indeed every solution to the constraints are physical.
Several comments are in order. First, similar constraints were obtained for gauge theories coupled to BSPT systems in Ref. 2 and Ref. 7. Most of the constraints here are just variants of those for BSPT phases, and some are even the same, e.g., Eqs. (8c), (8d) and (8g). However, Eq. (8e) is more "fermionic" than others, since it has no bosonic analog. It replaces the stronger condition N i Θ i,k = 0 in bosonic theories(see Ref. 7). Nevertheless, in a sense it is a "2D" constraint 25 , since the base loop does not enter the constraint in any nontrivial way. Similarly, Eq. (8i) has no analog in bosonic systems.
Second, the constraint Eq. (8h) remains a conjecture at this stage. We are not able to give a general proof. A weaker constraint can be derived from Eqs. (8a), (8f) and (8g): This weaker result provides some evidence for the conjecture Eq. (8h). In fact, it remains a conjecture in gauged BSPT systems too 7 (though, see Ref. 31 for a derivation of Eq. (8h) by exploiting the bulk-boundary correspondence under certain assumptions). Third, Eq. (8e) holds only for even N i with i ≥ 1. There are no analogous constraints for odd N i . In addition, Eq. (8f) holds only when N ik is even.
Lastly, we derive several useful corollaries. The first one follows from (8a) and (8c): Note that for BSPT theories, we have a stronger condition N ik Θ i,k = 0 (see Ref. 7). Another corollary follows from (8d) and (8e). Setting k = i in Eq. (8e) and using Eq. (8d), we immediately obtain Finally, combining Eq. (11) with Eq. (8f), we have Even though Eq. (12) follows from Eq. (11) which holds only for even N i with i ≥ 1, one can easily check that it holds in general.

C. Solutions from BSPT phases
Later on, we will solve the constraints Eqs. (8a)-(8i), where each solution leads to a consistent set of loop braiding statistics and corresponds to an FSPT phase. BSPT phases form a subset of FSPT phases, so we first write down a class of solutions that stem from BSPT phases.
Physically, we can imagine the following construction: first, we form "molecules" from pairs of fermions, where each molecule is a boson. The bosonic molecules are neutral under Z f 2 , and thus only sense the quotient sym- We then put the molecules into a BSPT state protected by the symmetry G. It is now generally believed that the 3D BSPT phases with unitary symmetry G are classified by the cohomology group H 4 [G, U(1)]. 32 The loop braiding statistics in gauged BSPT models were studied in Ref. 2 and 7. Now we can adapt the loop braiding statistics of BSPT phases from Ref. 7 into Θ i,k and Θ ij,k and find the following expressions: and Here M ijk is an arbitrary three-index integer tensor, N ik = gcd(N i ,N k ), and One can easily check that the above expressions satisfy all the constraints (8a)-(8i). Note that different values of M ijk can lead to the same values of the topological invariants.

D. Braiding statistics of fermion parity loops
We are mainly interested in loop braiding statistics beyond those given by Eqs. (13) and (14), i.e., those resulting from gauging "instrinsic" FSPT phases. Such loop braiding statistics will be explicitly discussed in the next subsection. Before doing that, we would like to answer this question: Given a solution to the constraints Eqs. (8a)-(8i), i.e., a set of three-loop braiding statistics, how do we know whether it is "intrinsically fermionic", and not just a gauged BSPT phase in disguise?
We claim that the following criterion holds: Criterion: A set of three-loop braiding statistics results from gauging an intrinsic FSPT phase, if and only if some of the three-loop braiding statistics involving fermion parity loops are "nontrivial".
Those three-loop braiding statistics that involve fermion parity loops include θ α,f , θ αβ,f , θ f ,γ and θ αf ,γ , where α, β, γ are arbitrary vortex loops and f stands for a fermion parity loop. We need to clarify what we mean by "nontrivial" threeloop braiding statistics. Just as any other vortex excitations, there are many fermion parity loops which differ by charge attachments. Attaching charges to loops shifts three-loop braiding statistics by Aharonov-Bohm phases. Accordingly, we call a three-loop braiding statistical phase "trivial" if it can be tuned to 0 by attaching charges to the loops involved in the braiding process. To remove the ambiguity due to charge attachment and Aharonov-Bohm phases, we define the following quantities for Abelian loop braiding statistics: where "lcm" stands for least common multiple, θ α,f is the exchange statistics of two identical α's linked to a fermion parity loop f , θ αβ,f is the mutual braiding statistics between α and β both linked to f , θ f ,e k is the exchange statistics of two identical f 's linked to a base loop γ, and θ f α,e k is the mutual statistics between α and f while both are linked to γ. Here, α, β, γ are vortex loops carrying gauge flux 2π Ni e i , 2π Nj e j , and 2π N k e k respectively. These quantities are defined in a way similar to the topological invariants Θ i,k and Θ ij,k . One can easily show that if Θ i,f , Θ ij,f , Θ f ,k and Θ f j,k vanish, all three-loop braiding statistics involving fermion parity loops are "trivial". Therefore to see if a set of three-loop braiding statistics corresponds to an intrinsic FSPT phase, we only need to check if any of the quantities Θ i,f , Θ ij,f , Θ f ,k and Θ f j,k is nonvanishing. This criterion can be proven by explicitly solving the constraints Eqs. (8a)-(8i) and checking if all the solutions with vanishing Θ i,f , Θ ij,f , Θ f ,k and Θ f j,k are in the form of Eqs. (13) and (14) (which will be discussed in the next subsection). On the other hand, Eqs. (13) and (14) indeed lead to vanishing Θ i,f , Θ ij,f , Θ f ,k and Θ f j,k .
Here, we would like to give a more intuitive argument. Although heuristic, it provides a physical interpretation for the criterion and can be applied more generally to non-Abelian loop braiding. Let us first show the "if" direction in the criterion. Consider those FSPT phases formed by bosonic pairs of fermions. Since the pairs do not transform under the fermion parity symmetry, they only need the symmetries in G/Z f 2 for protection. Under the assumption that FSPT phases have a one-to-one correspondence to equivalence classes of three-loop braiding statistics (i.e. up to charge attachment), it is reasonable to expect that the fermion parity loops should not play a nontrivial role in three-loop braiding statistics beyond Aharonov-Bohm phases after gauging the symmetry G. Hence, the "if" direction holds.
To see the "only if" direction, let us assume that fermion parity loops do not play any role beyond Aharonov-Bohm phases. Then, for any three-loop structure involving a fermion parity f , we can always attach charges to the loops such that (1) the new fermion parity loop f ′ has a bosonic exchange statistics and (2) for any given gauge flux, there always exists a vortex such that its mutual braiding statistics with respect to f ′ is trivial, while both are linked to the same base loop. Accordingly, we can condense f ′ , confining the fermionic gauge charges. This condensation leaves behind a gauge theory with purely bosonic charges, without affecting the other gauge symmetries, i.e., the resulting theory has a gauge group G/Z f 2 . As a consequence, the corresponding FSPT phases are always equivalent to those formed by bosonic pairs. This argument is reasonable but not quite rigorous, because a complete theory of loop condensation in topological phases is not available yet. We notice that similar arguments have been applied in two dimensions 33 .
Let us now combine this criterion with the constraints Eqs. (8a)-(8i) of topological invariants. We will see that the quantities Θ i,f , Θ ij,f , Θ f ,k are forced to vanish due to the constraints; only Θ f i,k may possibly be nonzero. Accordingly, for Abelian loop braiding statistics, we only need to compute Θ f i,k to see if a given set of Abelian three-loop braiding statistics corresponds to an intrinsic FSPT phase.
To derive these results, we first relate the quantities de-fined in Eq. (16) to the topological invariants as follows: These relations follow straightforwardly from the definitions of related quantities. The fact that Θ i,f = 0 follows immediately from Eq. (12), and Θ f ,k = 0 follows immediately form Eq. (8i). To see Θ ij,f = 0, we consider three cases.
Without loss of generality, we assume r i ≤ r j . It is not hard to see that Θ ij,f can only be 0 or π. Using Eqs. (8c) and (17b), we find that Θ ij,f may be nonzero only if r 0 ≤ r i ≤ r j . Assuming this is the case, we write Eq. (8g) in the following form: With this equation, we then have where the first line uses the facts that Θ ij,f can only be 0 or π and that t 0i t 0j t 0ij t ij is odd, the second line uses Eq. (19), the third line uses the constraint Eq. (8e), and the last line is a simplification. It is easy to see that both N i and N j divide the coefficient t 0i t 0j t 0ij 2 rj . Then, using the constraint Eq. (8d), we prove that Θ ij,f = 0.
Therefore, only Θ f i,k is potentially nonzero. Accordingly to Eqs. (8c) and (17d), one can see that it can only take values 0 or π. Moreover, if either N i or N k is odd, Eq. (8c) is enough to guarantee Θ f i,k = 0. If both N i and N k are even, from Eq. (8e), we have Θ f i,k = N i Θ i,k for i ≥ 1, which further lead to Θ f i,i = 0 using Eq. (8h). For i = 0, we have Θ f 0,k = N0 2 Θ 00,k = 0 according to Eq. (17c).
To summarize, we have shown that to check whether a set of topological invariants corresponds to an intrinsic FSPT phase or not, we only need to check if Θ f i,k is nonzero for i = k, i ≥ 1, and N i and N k are both even.

E. Solving the constraints
We now explicitly solve the constraints Eqs. (8a)-(8i). Mathematically speaking, the constraints are linear equations of the tensors Θ i,k and Θ ij,k . Solving them is straightforward, though tedious due to the fact that the equations are defined modulo 2π.
We first notice the following structure of the solutions: Given two sets of topological invariants Θ (1) and Θ (2) (2) . Due to the linearity of the constraints, Θ ′ also satisfy all constraints. In fact, we see that Combined with N i Θ ′ i,k = 0 we obtain N ik Θ ′ i,k = 0. In fact, Θ ′ satisfies essentially the stronger constraints for BSPT phases, whose solutions have been given in Sec. III C. Therefore, once we know the solutions corresponding to intrinsic FSPT phases, all others can be obtained by adding BSPT solutions.
To solve the constraints, a useful observation is that the constraints only relate those components of tensors whose indices differ at most by one index 0. Accordingly, we can divide the components of the tensors into four groups: where i = j = k = 0. Since Θ ij,k is symmetric in the first two indices, we do not list other components that are related by this symmetry above.
In the group (a), only the trivial solution is allowed: Θ 00,0 = Θ 0,0 = 0. It follows directly from the constraints Eq. (8a) and (8h). Also, invariants in the group (d) satisfy the same equations with those of BSPT phases. Hence, solutions for the group (d) are exactly the same as in BSPT phases, i.e., all can be written in the form Eqs. (13).
Below we solve the constraints for cases (b) and (c). Without loss of generality, we consider Consider the symmetry group G = Z f N0 × Z N1 with N 0 = 2m. We solve the constraints for topological invariants in the group (b) with i = 1. Among the eight components, we find that Θ 0,1 , Θ 1,0 and Θ 1,1 completely determine the rest. More explicitly, where the first and second lines follow Eq. (8a) and the third line follows Eq. (8f). Since Θ 1,1 = 0 following Eq. (8h), we only need to find possible values for Θ 0,1 and Θ 1,0 . Let us first consider odd N 1 . In this case, using the constraints Eqs. (8a),(8c), (8d), and (8i), we find that Here x, y are integers. It is not hard to see that this solution is in the form of Eq. (14). This agrees with the criterion discussed in Sec. III D.
When N 1 is even, using Eqs. (8d) and (8i), we find that Θ 0,1 can take the following values where x is an integer. According to the corollary Eq. (12), we have mΘ 1,0 = 0. In addition, multiplying Eq. (8e) by 2 and using Eq. (8c), we have 2N 1 Θ 1,0 = 0. Together we find For odd m, the parameters x and y are independent. For even m, there exists a relation between x and y: Taking i = 1 and k = 0 in Eq. (8e) and using the expression of Θ 01,0 = Θ 10,0 in Eq. (22), we find that This relation puts an constraint x ≡ y (mod 2) on x and y only if N 1 / gcd(m/2, N 1 ) is odd.
Let us see which solutions correspond to "intrinsic" FSPT phases. According to the criterion discussed in Sec. III D, we only need to check the quantity Θ f 1,0 . It is non-vanishing only when m is even and N 1 is even, in which case we find that Let N 0 = 2 r0 t 0 and N 1 = 2 r1 t 1 , where r 0 ≥ 2, r 1 ≥ 1 and t 0 , t 1 are odd numbers. Then, it is easy to see that Θ f 1,0 = π only if and x is an odd number. Therefore, the simplest symmetry group to support intrinsic FSPT phases is In this case, one will find that if either N 1 or N 2 is odd, all solutions to the constraints correspond to BSPT phases, given by (13) and (14). Hence, below we focus on the more interesting case where both N 1 and N 2 are even. One can show that any odd factors of N 0 , N 1 , N 2 cannot add solutions that correspond to intrinsic FSPT phases. Hence we assume that N i = 2 ri for i = 0, 1, 2, with r 0 , r 1 , r 2 ≥ 1 for simplicity. Without loss of generality, we further take r 1 ≤ r 2 .

F. Examples
In this subsection, we discuss two examples whose three-loop braiding statistics correspond to intrinsic FSPT phases.
One of the simplest groups that support intrinsic FSPT phases is G = Z f 8 × Z 2 . It supports the following threeloop braiding statistics and all other invariants are 0. It is a solution obtained in Sec. III E by taking x = y = 1 in Eqs. (24) and (25).

Physical picture
Let us understand the two examples in more physical terms. Although the two examples have seemingly different symmetry groups, they are in fact closely related. Both symmetry groups can be regarded as central extensions of Z 2 × Z 4 : namely, we can take the physical symmetries to be G = Z 2 × Z 4 in both cases, and for Z f 8 × Z 2 the fermions carry half charges under the Z 4 subgroup.
Consider a 2π 4 e 2 base loop. By dimensional reduction we obtain a 2D fermionic SPT protected by the Z 2 × Z 4 symmetry. In fact, because of Θ 1,2 = π 2 , the protecting symmetry is just the Z 2 subgroup. As we have already discussed in Sec. I, this is the "root" Abelian Z 2 FSPT phase in 2D, which has a Z 4 classification (thus can only exist on a Z 4 base loop) 14,20,22 . Besides the noninteracting realization mentioned in Sec. I, commutingprojector Hamiltonians for such 2D phases have also been found in Ref. 14, 34, and 35. Θ 1,2 = π 2 translates into the fractional exchange statistics of symmetry fluxes in the 2D FSPT phase, which as proven in Ref. 20 implies its edge (i.e. the Z 4 vortex loop) has to be degenerate.

IV. TOPOLOGICAL STATE-SUM MODELS
In this section we introduce a class of lattice models to realize the fermionic gauge theories found in the previous section. We define these lattice models with a path integral representation of the partition function in discretized Euclidean space-time. More specifically, we define a partition function for any closed oriented manifolds with a triangulation. The partition function however is a topological invariant of the space-time manifold (i.e. independent of the choice of the triangulation). Hence, it is a type of lattice topological quantum field theory (TQFT). It is generally believed that such topological state-sum models can be cast into commuting-projector Hamiltonians 36 .
We will first recall a few useful facts regarding triangulations of n-dimensional manifolds. We will work with simplicial triangulations for simplicity 37 and denote the set of k-simplices (0 ≤ k ≤ n) in the triangulation as ∆ k . For a given triangulation, we pick an arbitrary ordering of the vertices as 0, 1, 2, . . . . The restriction of the ordering on each k-simplex σ k induces a relative ordering of the vertices of σ k . Under this relative ordering, we write σ k as [i 0 i 1 · · · i k ], where i 0 < i 1 < · · · < i k are the vertices of σ k .
On an oriented manifold, all simplices can be equipped with orientations, induced from the orientation of the manifold M . For each σ n , we define ε(σ) to be 1 if the orientation on σ induced from that of M coincides with the one coming from the ordering of its vertices; otherwise if they are opposite then ε(σ) = −1.

A. Twisted Crane-Yetter TQFT
We now present models for fermionic gauge theories. The construction was first introduced by Kapustin and Thorngren recently in the context of higher-form gauge theories 38 . We will call these models the twisted Crane-Yetter models. The general input of the twisted Crane-Yetter TQFT involves (i) a braided fusion category (BFC), (ii) a finite group G and (iii) certain cohomological data (β, λ, ω) associated with G and the BFC.
For simplicity, we present the construction for an Abelian BFC A. The anyon labels in A are denoted by a, b, c, . . . . The identity (i.e., the trivial anyon) is denoted by 0. The BFC A can be viewed as an Abelian group, with the group multiplication given by the fusion rules. As a BFC, A is equipped with further topological data, in particular the F and R symbols. We will further assume that the F symbols of A can be chosen to be trivial. In this case, the hexagon equations simplify to Notice that because A is Abelian, we denote the multiplication additively. In other words, R a,b defines a bicharacter on the Abelian group A. We define T as the following subset of A: We refer to T as the subset of transparent particles. In many cases, we will actually take a BCF such that A = T . The other pieces of the input data are a finite group G and two group cocycles: a 3-cocycle [β] ∈ H 3 [G, T ], and a 2-cocycle [λ] ∈ H 2 [G, A] ([·] denotes the cohomology class). A 3-cocycle β is a 3-cochain (i.e., a function) G 3 → T that satisfies the 3-cocycle condition: Similarly, a 2-cocycle λ is a function G 2 → A that satisfies the 2-cocycle condition: The final piece of data ω is a group 4-cochain ω : G 4 → U(1). It is generally speaking not a 4-cocycle, however it does satisfy a similar condition which will be determined later. We will frequently use the following short-hand notation for a group n-cochain ν: face, i.e. the connection is flat. So for each 2-simplex [ijk] (i < j < k) the flatness condition is imposed: To each 2-simplex [ijk] (i < j < k), we assign a simple object f ijk from A. For each 3-simplex [ijkl], we demand that the following "flatness condition" holds: Let us now write down the partition function. To each 4-simplex, say σ 4 = (01234), we assign a phase factor: The phase factor T + (σ 4 ) is assigned to σ 4 assuming its local orientation coincides with the global orientation of M . If they have the opposite orientations, we instead assign T − (σ 4 ) = T + (σ 4 ) * to σ 4 . The partition function is then defined as We require that Z(M ) defines a topological quantum field theory. Namely, Z(M ) should yield a topological invariant of the manifold M , which means that 1) it must be independent of the specific choice of triangulation and 2) independent of the ordering of the vertices. It is known that all triangulations can be related to each other via a finite series of elementary moves, known as Pachner moves 39 . To this end, we define the following "obstruction class" One can show that O is actually a 5-cocycle in H 5 [G, U(1)]. We show in Appendix D that invariance under Pachner moves requires that We observe that the right-hand side of Eq. (47) is the coboundary of the 4-cochain ω. Hence, it implies that the obstruction class must be cohomologically trivial in order for the twisted Crane-Yetter model to be well defined. Otherwise, we say the model is "obstructed". For obstruction-free models, Eq. (47) becomes a "twisted" 4cocycle condition on ω (compared to the regular 4-cocycle condition in which the left-hand side is 1). The twisted Crane-Yetter models reduce to known models in two special limits: 1. G is trivial. In this case, the state sum reduces to the well-known Crane-Yetter theories 40,41 (the Hamiltonian version of the TQFT is known as the Walker-Wang model 42 in the condensed matter literature). Excitations in the model can be understood as a T gauge theory, however with an interesting twist: particle excitations are labeled by elements of T . A particle a then has topological spin θ a = R a,a = ±1. In fact, the characterization of particle excitations hold generally, not just for the Abelian BFCs discussed here. Therefore, in general Crane-Yetter models also produce topological gauge theories. Recently it is shown that with non-Abelian BFC as the input the Crane-Yetter model can also realize twisted gauge theory. 43

2.
A is trivial. In this case, the theory reduces to the Dijkgraaf-Witten topological gauge theory 44 .

B. Gauge-theoretical interpretation
The state-sum model can be understood as a topological gauge theory, for a 2-form gauge field f and a 1-form gauge field g. In other words, the theory embodies two kinds of gauge symmetries: the 1-form gauge transformations on f and 0-form gauge transformations on g: We explicitly prove that the partition function is invariant under the two types of gauge transformations, with two simplifying assumptions: (a) G is Abelian and (b) The entire A is transparent (i.e. T = A). While the proof is rather technical and the details can be found in Appendix E, we introduce the reformulation of the TQFT as a topological gauge theory following Ref. 38, using the notations of simplicial calculus (see Appendix C for a review of relevant mathematical concepts). In the following the multiplication in G will also be denoted additively. We define the discretized "action" S by T(σ 4 ) = e iS(σ4) . f can be viewed as a 2-cochain valued in A, and g is a 1-cochain valued in G. The flatness conditions Eqs. (42) and (43) then can be written as δg = 0 and δf = β. The latter implies that f is not closed.
In the partition function, the product of three R symbols closely resembles the "Pontryagin square" in Ref. 38.
Roughly speaking, if f is a closed 2-cochain, the Pontryagin square is just the cup product f ∪ f . However, if f is not quite closed, the cup product is no longer closed and we need to amend it by an additional term: f ∪f −f ∪ 1 δf to get a closed cochain. In this notation, the action can be written as Here η = ln ω 2πi is the linearized form of ω. With this notation, we can now give a quick derivation of the obstruction-vanishing condition. In order for the action S to be a topological invariant, all we need to show is that S is a closed 4-cochain: We thus require δη = β ∪ 1 β + λ ∪ β, which is the obstruction-vanishing condition Eq. (47).

C. Relation to symmetry-enriched topological phases
We now define a variant of the Crane-Yetter TQFT. Instead of having G labels on the 1-simplices, we dualize them to 0-simplices, i.e. vertices. Namely, the actual labels are G group elements on vertices, and g ij = g −1 i g j . The flatness condition for f 's is the same as before, as well as the partition functionT + (σ 4 ): Here T + (σ 4 ) is the partition function defined in Eq. (44), where the G label on [ij] is given by g −1 i g j . This "dual" state-sum model can be viewed as a model of symmetry-enriched topological phase. For each h ∈ G, a global symmetry transformation is defined as Apparently the partition function is invariant under such global symmetry transformations since it only depends on g −1 i g j . As we have argued, the G fields can be understood as connections of a discrete G gauge field. As a result, gauge-equivalent configurations of G fields yield the same partition function. In the SET model, the connections {g ij } are by definition "pure gauge". Thus the partition function on any oriented closed manifold is identical to g i = 1, i.e. the Crane-Yetter TQFT of A. The original state-sum model can then be viewed as "gauging" the SET model. This relation was first considered between DW gauge theories and group-cohomology models of bosonic SPT phases 18 .
Therefore, we can understand that the bulk excitations of the SET model are described by the Crane-Yetter TQFT, i.e. a T gauge theory (with possibly fermionic gauge charges). An important question is then how excitations, including both particles and loops, transform under the symmetry group G. In the next section, using dimensional reduction, we will argue that the particles transform as projective representations of G essentially determined by λ, while the loop excitations exhibit nontrivial symmetry actions corresponding to nontrivial three-loop braiding statistics after gauging.

V. BRAIDING STATISTICS IN THE STATE-SUM MODELS
In this section, we analyze the braiding statistics between particle and loop excitations in the twisted Crane-Yetter state-sum models. For simplicity, we assume that the whole BFC A is transparent, i.e., A = T , throughout the section.

A. Particle excitations
First let us consider the properties of particle excitations. For this purpose, it is useful to view the state-sum model as a G-symmetry enriched gauge theory. As we have argued earlier, the particle excitations are charged under the emergent gauge group A. We need to understand how they transform under G. We will present some evidence that the symmetry transformations on particles are determined by λ.
To begin, we consider a simpler theory defined by together with the twisted flatness condition δf = β. Basically we drop the term in the action which contributes to the exchange statistics of the particle excitations (so now they are all bosons). Ref. 38 showed that this theory can be analyzed in a dual representation, where the 2-form is dualized to a 1-form gauge field a. To simplify the discussion, we use the simplicial calculus. First the constraint δf − β = 0 can be implemented by introducing a 1-cochain Lagrange multiplier a and the following modification to the action: We should notice that this applies only to non-degenerate R, namely 1 |A| a∈A R a,b = δ(b). Since f is no longer constrained, we can sum over f to get δa − λ = 0, and the action becomes S = −a ∪ β. This is a symmetryenriched A gauge theory, where the gauge charges, labeled by characters χ of A, transform as projective representations of G. More specifically, the projective representation U χ on the gauge charge χ is given by where the projective phase factor is given by In general, with the f ∪ f term in T + (σ 4 ), we can not apply the above duality transformation. However for a special case A = Z f 2 and with a trivial β, we can "linearize" the term f ∪ f using the following relation: where w 2 is the second Stiefel-Whitney class of the manifold. This relation only holds when f is a 2-cocycle. Then a similar duality transformation leads to δa = w 2 + λ. The w 2 term accounts for the fermionic statistics of particles 45 , and the symmetry transformation under G is again given by Eqs. (56) and (57). Extrapolating from these two special cases, we conjecture that Eqs. (56) and (57) hold more generally in the full theory (when the R symbol is nondegenerate). We do not have a proof of this statement at the moment, but will show that it is also consistent with the dimensional reduction result.

B. Dimensional reduction
To understand properties of the loop excitations, we consider the theory on a 4-manifold M 4 = M 3 × S 1 , where M 3 is a 3-manifold and S 1 is the circle. Since these models have zero correlation length, we can analyze the theory in the limit where there is only one cell in the S 1 direction, and view it as a 2D theory on M 3 . We will fix the G flux through the "hole" of S 1 to be h.
Some of the particle excitations in the 2D theory correspond to the loop excitations that are linked to a h gauge flux in the 3D theory. Accordingly, if we can extract the braiding statistics in the 2D theory, three-loop braiding statistics in the 3D theory can be inferred from there.
Let us understand the fields in the 2D theory after dimensional reduction. All the 1-form G gauge fields in M 3 remain, so do the 2-form A gauge fields in M 3 . Both of them satisfy the same (twisted) flatness conditions as before. However, there are additional dynamical 1-form fields, denoted by m, coming from the dimensional reduction of the full 2-form gauge fields in M 4 . They obey the following flatness conditions: Here n is the slant product of β: n = i h β. The explicit expression of n in terms of β reads: The partition function with a fixed holonomy h around S 1 will be denoted by Z h (M 3 × S 1 ). After a straightforward but lengthy calculation, we find In this expression ∆ k ≡ ∆ k (M 3 ), and the weight on a tetrahedron is given by Here the 1-cocycle ξ is the slant product of λ: and α is a 3-cochain that depends on β and ω. The obstruction vanishing condition Eq.(47) reduces to the following equation We refer the readers to Appendix F for the derivation of the dimensional reduction.

C. Braiding statistics in the 2D theory
After dimensional reduction, we now analyze the 2D theory. To simply the analysis, we assume that G is Abelian, and ξ in Eq. (63) is cohomologically trivial. With the latter assumption, we see that in the 2D partition function Z h the sum over f becomes completely independent of h, and the h dependence only enters through the sum over m. Therefore for the purpose of extracting loop braiding statistics on the h base loop, we only need to focus on the sum over m. More discussions on the physical meaning of the sum over f can be found in Appendix F. It turns out that to realize those threeloop braiding statistics studied in Sec. III, it is sufficient to assume ξ is cohomologically trivial.
To infer the loop braiding statistics in the original 3D theory, our analysis of the 2D theory should proceed in two steps. First, we need to establish the correspondence between the excitations in the 2D and 3D theories. More precisely, we need to identify which of the 2D anyons correspond to the 3D particle excitations and which correspond to the 3D loop excitations. Second, we need to extract the braiding statistics of the 2D anyons.
All properties of the 2D theory Z h should depend on h. For notational brevity, we suppress this h dependence below. It is easy to recover this dependence later.

Correspondence between 2D and 3D excitations
As discussed above, the 2D theory has two kinds of dynamical variables, m ij ∈ A and g ij ∈ G living on each link [ij]. On each 2-simplex [ijk], they satisfy the twisted and untwisted flatness conditions respectively: Instead of viewing m ij and g ij as two independent degrees of freedom, we can combine them into one and denote it as (m g ) ij . In fact, {m g }| m∈A,g∈G form a group G under the following definition of group multiplication: The groupG is known as a central extension of G by A, associated with the 2-cocycle n(g, k) ∈ H 2 [G, A].
With this notation, we claim that the partition function T 2+1 actually represents a 2D Dijkgraaf-Witten gauge theory of groupG associated with the following 3-cocycle To see this, one can check that (i) ω 2+1 is indeed a 3cocycle in H 3 [G, U(1)], as long as n, λ, α satisfy Eq. (64) and (ii) the conditions (65) leads to i.e., every 2-simplex has a flatG connection. The above mapping is most convenient for general computations of braiding statistics of excitations in the 2D theory, since braiding statistics in Dijkgraaf-Witten theory is known (see e.g. Ref. 46). Below we will take a slightly different approach. We will dualize the G gauge fields in theG Dijkgraaf-Witten theory, similarly as discussed in Sec. IV C, and view the result as an A gauge theory, enriched by the symmetry group G.
The Hamiltonian version of this symmetry-enriched gauge theory was recently considered in Ref. 47. The anyons in an Abelian A gauge theory can be labeled as dyons (a, χ) where a ∈ A is the gauge flux, and χ : A → U(1) denotes a character of A, labeling a gauge charge. Since the symmetry group G does not permute any anyons, each anyon carries a projective representation of G. The twisted flatness condition Eq. (65) is interpreted as A gauge charges carrying projective representations of G. As shown in Ref. 47, the projective phases on a pure charge (1, χ) is For a gauge flux (a, 1), η (a,1) (g, k) = R a,n(g,k) R λ(g,k),a More generally, the projective phase of the dyon (a, χ) is given by η (a,χ) = η (a,1) η (1,χ) . We now identify the correspondence between the 2D and 3D excitations. For each a, we define χ a (x) = R x,a with x ∈ A. It is clear that χ a is a character of A. Then, the dyon (a, χ a ) transforms under G with a projective phase η (a,χa) (g, k) = R a,n(g,k) R λ(g,k),a R n(g,k),a = R λ(g,k),a = χ a λ(g, k) (71) where we used the fact that β, and therefore n, is transparent to obtain the second line. It is easy to see that the anyons {(a, χ a )} a∈A form a fusion group A, and we identify them as the descendants of the 3D quasiparticles. Indeed, the topological twist of (a, χ a ) is χ a (a) = R a,a = θ a , and the mutual braiding between (a, χ a ) and (b, In fact, it is not difficult to see that the set {(a, χ a )} forms the maximal subset of transparent anyons in the 2D theory, i.e. any other particles not in this set have nontrivial braiding with at least one of the anyons in the set. As a result, all other particles should be understood as the descendants of gauge flux loops in the 3D theory.

Braiding statistics with
With the above understanding, we now specialize to the case A = Z N0 and G = K i=1 Z Ni . For the rest of the section we assume the even number N 0 = 2m, and the R symbols are given by R a,b = (−1) ab , which corresponds to Z f N0 . We will also work with the followng parametrization of λ: where q i are integers, and [g i + k i ] equals g i + k i modulo N i . We have used integer vector g = (g 1 , . . . , g K ) to deonote group elements of G. In fact, this parameterization exhausts all cohomology classes in H 2 [G, Z N0 ] which satisfy the assumption that ξ in Eq. (63) is cohomologically trivial. We would like to extract a part of the braiding statistics data from given n(g, k) and α(g, k, l), focusing on those that will lead to Abelian statistics. To start, let us give an explicit parametrization of the inputs n(g, k) and α(g, k, l). Let us take the following class of 2-cocycles n ∈ H 2 [G, Z N0 ]: where p i are integer parameters. Here, we have used additive convention for group multiplication in both G and A. It is worth to mention that this class of n satisfies n(g, k) = n(k, g). This property is a necessary and sufficient condition for the braiding statistics to be Abelian. At the same time, we can choose where t ij are integers. One can check that the above λ, n, α indeed make ω 2+1 a 3-cocycle. The fact the existence of α for given λ and n (not just those parametrized by Eq. (73)) shows that there is no obstruction in the 2D theory. Notice that if we modify the 4-cochain ω by a 4-cocycle, α will be modified by a 3-cocycle correspondingly.
One can in principle compute the braiding statistics using the understanding from previous subsection and using the general results for Dijkgraaf-Witten models 46,48 , but in our case we will take a short cut. Notice that in our parametrization of n(g, k), different Z Ni subgroups are decoupled, so let us focus on an individual Z Ni for now. In the 2D Z N0 gauge theory, following the discussion of dual SET phases in the previous subsection we denote the unit gauge flux by v, and the unit gauge charge e corresponding to the character λ e (x) = exp( 2πi N0 x), with x ∈ A. The (bulk) fermion in this notation is represented by e N0/2 v.
With a Z Ni global symmetry, anyons can carry fractional charges under Z Ni . Denote the projective phase on an anyon a by η a (g, k). We can extract the fractional charge Q a as follows: .
We can say that e anyon carries a pi N0i fractional charge of the Z Ni symmetry, while v anyon carries N0(pi+qi) 2N0i fractional charge. The fermion e N0/2 v carries a fractional charge N0qi 2N0i , which can only be 0 or 1 2 . Such an SET can be easily described by an Abelian Chern-Simon theory with the following K matrix and charge vector t i 27,49,50 : In this formalism, anyons are labelled by integer vectors. The exchange statistics θ l of an anyon l and mutual braiding statistics between anyons l and l ′ are given by In addition, the Z Ni charge carried by l is given by We denote the gauge charge by e = (1, 0) and the gauge flux v = (0, 1). One can easily check that such K and t i indeed descirbe the above Z N0 gauge theory enriched by Z Ni symmetry. So far we have focused on the symmetry fractionalization in the 2D theory, which is completely determined by β and λ. We have not accounted for the possibility of adding 2D BSPT layers (i.e. different 3-cocycle α). While the details of braiding statistics surely depend on the BSPT layer, this subtlety does not affect the indicator Θ f i,k for intrinsic FSPT phases, discussed in Sec. III D, which is of central interest to us. So, we ignore the BSPT ambiguity (i.e., α dependence) for the moment.
We now gauge the symmetry G and switch back to the twisted Crane-Yetter model. After gauging, a vortex carrying the unit Z Ni flux can be represented by the fractional vectors t i /N i . The braiding statistics between these vortices can be calculated again using Eq. (78). Accordingly, we obtain the exchange and mutual braiding phases: The mutual braiding between the Z Ni unit flux and the e anyon, as well as the v anyon, is given by These braiding statistics are only part of the full set of braiding statistics. Note that there exist other Z Ni vortices differing by charge attachments. However, knowing the braiding statistics between the vortices {t i /N i } i=1,...,K , as well as the anyons e and v, is enough to extract the topological invariants for gauged intrinsic FSPT phases defined in Sec. (III).

D. Braiding statistics in the 3D theory
Now we come back to the 3D theory. Again we consider G = K i=1 Z Ni and A = Z f N0 , with R a,b = (−1) ab . The 2-cocycle λ is given in Eq. (72). We will ignore the ω dependence of the braiding statistics, which is not relevant for intrinsic FSPT indicator Θ f i,k .
We use the following explicit expression for representative cocyclces β ∈ H 3 [G, Z N0 ]: where p ij and p ijk are integer parameters, and [h j + k j ] equals h j + k j modulo N j . It turns out that the first term in the explicit expression describes Abelian braiding statistics for loops, while the second term leads to non-Abelian loop statistics. For this reason we will refer to the two types of cocycles as "Abelian" and "non-Abelian". In this subsection, we only consider the Abelian part of the cocycle, leaving a discussion on non-Abelian loop braiding statistics in Sec. VII. One may calculate the cohomology group H 3 [G, Z N0 ] using the Künneth decomposition: We believe that the parametrization in Eq. (82) exhausts all cohomology classes in H 3 [G, Z N0 ]. We now focus on the Abelian part of β. As discussed previously, dimension reduction of the 3D statesum model leads to a 2D model described by the slant product n = i h β. For β in Eq. (82), the slant product is given by where This form of n is the same as in Eq. (73), with p i there replaced by P i . Taking h = e k and substituting the expression of P i into Eqs. (80) and (81), we obtain the three-loop braiding statistics: and Previously, we have identified e N0/2 v as the fermion particle in the 3D bulk. Since e N0/2 and v both have a π mutual statistics with respect to e N0/2 v, we now should understand both of them as fermion parity loops. Accordingly, we understand v 2 as a bosonic particle in the 3D bulk, and e is a loop excitation that has a 2π N0 mutual statistics with respect to the fermion particle e N0/2 v.
Finally, we make the following comment. Throughout our computation, we do not keep track of the dependence of three-loop braiding statistics on the 4-cochain ω. A consequence is that the three-loop braiding statistics given in Eqs. (86) and (87) are not the complete result, in particular violating the constraints in Eqs. (8e) and (8f). In order for the constraints to be obeyed, we have to keep track of the ω dependence carefully, which however is very complicated. Nevertheless, the indicator Θ ℧i,k of intrinsic FSPT phases does not change after attaching BSPT layers, therefore we can safely ignore the issue for the purpose of extracting the indicators.

E. Realizations of intrinsic FSPT phases
We now show that the state-sum model realizes all intrinsic FSPT phases that we found in Sec. III, completing the argument that the physical constraints Eqs. (8a)-(8i) are complete. In Sec. III, we find two kinds of intrinsic FSPT phases, supported by the representative groups Z f 2m × Z N1 and Z f 2m × Z N1 × Z N2 respectively. Without loss of generality, we assume m = 2 rm , N 1 = 2 r1 and N 2 = 2 r2 . According to Sec. III E, existence of the two kinds of intrinsic FSPT phases requires r m ≥ r 1 + 1 ≥ 2 and r 2 ≥ r 1 + 1 ≥ 2 respectively. Since the second kind of intrinsic FSPT does not put requirements on m, we assume m = 1 for simplicity. One can easily extend the following discussion to general m for the second kind of intrinsic FSPT phases.
Let us now take a unified view on the groups Z f 2m ×Z N1 and Z f 2 × Z N1 × Z N2 : both of them arise as central ex- More specifically, the former is a nontrivial central extension of Z N1 × Z N2 by Z f 2 , associated with N 2 = m and q 2 = 1, q 1 = 0 in the 2-cocycle λ, while the latter is the trivial extension Z N1 × Z N2 by Z f 2 associated with q 2 = q 1 = 0 in λ. Therefore, in the twisted Crane-Yetter state-sum model, we set A = Z f 2 and choose λ accordingly. 51 We will see below that through this choice, the state-sum model indeed can be viewed as Z f 2m × Z N1 and Z f 2 × Z N1 × Z N2 gauge theories coupled to fermionic matter respectively. In this notation, the condition r m ≥ r 1 + 1 ≥ 2 on existence of intrinsic FSPT phases for Z f 2m × Z N1 translates to r 2 ≥ r 1 + 1 ≥ 2, the same as that for Z f 2 × Z N1 × Z N2 . Let us specify the input data to the state-sum model. As discussed above, for the 2-cocycle λ in Eq. (72), we set q 1 = 0, and q 2 = 0, 1 for Z f 2 × Z N1 × Z N2 and Z f 2m × Z N1 respectively. For both groups, we set the 3-cocycle β in Eq. (82) by the following parameters: The non-Abelian part of β is set to 0. Before we discuss the loop braiding statistics, we need to check that with these choices of λ and β the topological state-sum model is well-defined, i.e. the obstruction class (46) vanishes. In Appendix G, we provide a complete set of invariants to distinguish all cohomology classes in H 5 [G, U(1)] when G is a finite Abelian group. Applying these invariants to the present case, we find that the obstruction class vanishes, when the following equations hold: where the second equation should hold for i = j, and N 0 = 2. These equations are defined modulo 2π. With our choice of the parameters q i , p i and p ij in λ and β, we find that the first equation in Eq. (89) is automatically satisfied, while the second one puts the following conditions on r 1 and r 2 : which does not depend on q 2 . For r 2 ≥ r 1 + 1 ≥ 2, we see that the above equation indeed holds modulo 2π. Hence, the twisted Crane-Yetter state-sum model is obstructionfree.
The loop braiding statistics for the two groups are given by Eqs. (86) and (87), under the current choice of q i , p i and p ij . Let us check that the braiding statistics imply that they are indeed Z f 2m ×Z N1 and Z f 2 ×Z N1 ×Z N2 gauge theories. Since N 0 = 2, the fermionic particle is ev, and the fermion parity loops are e and v. According to Eq. (87), the mutual braiding between ev and the Z N2 unit flux on any base loop is given by q 2 π/N 2 . The mutual braiding statistics between ev and the Z N1 unit flux is always 0. Hence, it is indeed a Z f 2m × Z N1 gauge theory for q 2 = 1, and a Z f 2 × Z N1 × Z N2 for q 2 = 0. With this understanding, we now calculate the indicator Θ fi,k for intrinsic FSPT phases: where we understand that e is the fermion parity loop, and N 0 = 2, N 1 = 2 r1 , N 2 = 2 r2 with r 2 ≥ r 1 + 1 ≥ 2. One may use the v fermion parity loop to do the computation, which leads to the same result. This agrees with the results in Sec. III E (the index "2" should be understood as "0" for the group Z f 2m × Z N1 ). Therefore, all intrinsic FSPT phases identified in Sec. III are realized in the twisted Crane-Yetter model.

VI. ANOMALOUS SETS IN 3D
In Sec. III we derived a set of physical constraints for Abelian loop braiding statistics. We now demonstrate that these constraints can be used to show that certain 3D symmetry-enriched gauge theories are anomalous.
Let us discuss a simple example: G = Z 2 and A = Z f N0 with N 0 = 2 and R a,b = (−1) ab . In a Z f 2 gauge theory enriched by Z 2 symmetry, if the fermion parity flux loop carries gapless modes whose symmetry transformations are identical to those of a nontrivial 2D Z 2 BSPT phase 18 , such a SET is anomalous. To see the anomaly, gauging the Z 2 symmetry we would obtain Θ 1,0 = π forbidden by the constraints. More specifically, Θ 1,0 = π implies Θ 01,1 = π, contradicting Eq. (8e).
However, such a SET can actually be realized consistently on the surface of a 4D bosonic Z 2 SPT phase, for instance, by a coupled-"layer" construction as presented in Ref. 52. Based on heuristic field theory arguments, Ref. 52 also proposed that the bulk 4D SPT state is the one obtained from group-cohomology classification (H 5 [Z 2 , U(1)] = Z 2 ). We can provide a more rigorous justification with the topological state-sum models. Let us set λ = 0 for the moment, and the topological action has a variation where ∼ means up to a 5-coboundary. In order for the model to be a well-defined topological gauge theory in 3D, δS has to vanish cohomologically. When the obstruction class does not vanish, we have to couple the model to a 4D theory. The fields in the bulk are just G spins on vertices, and the topological action is given precisely by β ∪ 1 β. Therefore the bulk is essentially a group-cohomology model of a bosonic SPT phase. Back to the example, let us take G = Z 2 , A = Z f 2 , and a nontrivial 3-cocycle given by β(g, g, g) = [1] (we represent A = {[0], [1]}). One can easily check that the obstruction class is nontrivial. If we naively apply the dimensional reduction method to compute loop braiding statistics, we would find Θ 01,1 = π.
We also notice that the same obstruction appears in the gauge theory with all bosonic charges Eq. (54) if we have λ(g, g) = [1] and the same β. For d = 3, the mathematical structure of the Gu-Wen construction is completely identical to the twisted Crane-Yetter TQFT with A = Z f 2 and a trivial λ. We believe that the state-sum model discussed in this work with λ = 0 is indeed a gauged Gu-Wen model, where fermions are coupled to Z 2 gauge fields. Our results also clarify the physical meaning of [β] ∈ H 3 [G, Z 2 ] for Abelian unitary G, that is, the cocycle β encodes information about the three-loop braiding statistics.
With a nontrivial λ the state-sum model generalizes the Gu-Wen supercohomology constructions, by allowing gauge fermions to carry projective representations of the symmetry group. We have considered "Abelian" cocycles for λ. It will be interesting to explore the physics of "non-Abelian" 2-cocycles, corresponding to fermions carrying higher-dimensional projective representation of the symmetry group. A recent discussion on such terms in continuum field theories can be found in Ref. 53.

B. Non-abelian loop braiding statistics
We have exclusively focused on Abelian loop braiding statistics in this work. Loops can also exhibit non-Abelian braiding statistics. This can happen even when the gauge group is Abelian, if we choose a "non-Abelian" 3-cocycle β in (82). We will present one such example, , λ = 0 and the 3-cocycle β is parametrized by p 123 = 1 with all other components of p set to 0. Using the invariants given in Appendix G, it is easy to show that the state-sum model is obstructionfree.
To see the non-Abelian loop braiding, consider a base loop φ 1 . From the dimensional reduction, the e and v anyons in the 2D theory both carry two-dimensional projective representations of G. After G is gauged, they become non-Abelian anyons and exhibit non-Abelian braiding statistics, similar to what has been found in certain Dijkgraaf-Witten gauge theories 7,8,54 .
Recent works have constructed exactly-solvable lattice models for putative non-Abelian 3D topological phases 55,56 , in which the dimensionally reduced theories may support non-Abelian Ising excitations. It will be interesting to extend the dimensional reduction approach to these models. We treat the topological invariants and three-loop braiding statistics interchangeably throughout the paper. Here, we show that they are indeed equivalent in the case of Abelian braiding statistics. The following argument is a simple generalization of that for BSPT phases given in Ref. 7.
To show the equivalence, it is enough to reconstruct the full set of three-loop braiding statistics out of the topological invariants. Consider an arbitrary set of vortices {v i } with v i carrying unit flux 2π Ni e i . All of them are linked to a base loop that carries unit flux 2π N k e k . According to the definitions of Θ ij,k and Θ i,k , we have where x ijk , y ik are some integers that satisfy the relations y iik = 2x ik and y ijk = y jik . These relations follow from the properties θ αα,γ = 2θ α,γ and θ αβ,γ = θ βα,γ . We take {Θ ij,k , Θ i,k } in the interval [0, 2π), but in certain cases we set some of them in the interval [2π, 4π), which will be discuss below. Then, we attach a charge q ik to the loop v i when it is linked to 2π N k e k unit flux, for each i and k. The new vortex loops {v i } have the following mutual and self three-loop braiding: where q ik j is the jth component of q ik . We choose the charge {q ik } properly such that they satisfy the following relations: One can show that for even N i , such {q ik } always exist. For odd N i (i ≥ 1), the existence of such {q ik } requires y 0ik = x ik (mod 2). Interestingly, it is actually a physical requirement for properly chosen Θ i,k and Θ i0,k . Before we explain the case of odd N i , we conclude that if Eq. (A3) holds, we obtain a set of vortex loops {v i } such that That is, these braiding statistics are determined by the topological invariants. We now explain the case when N i is odd for i ≥ 1. Consider N i copies of v i vortices linked to 2π N k e k flux. Fusing the v i vortices together gives a pure charge q. Using the linearity relations (A8c) and (A8d) which we will explain shortly, we obtain (More detailed discussion can be found in Appendix B in the proof of Eq. (8e).) In addition, using the constraints Eqs. (8a), (8c) and (8d) for odd N i , we can write the topological invariants as follows

Inserting Eqs. (A1) and (A6) into
Eq. (A5), we find that π(a ik + x ik ) = π(b ik + y 0ik ), (mod 2π) (A7) Now if a ik and b ik have the same parity, so do x ik and y 0ik . If a ik and b ik have opposite parity, we can replace a ik by a ik + N ik in Eq. (A6). This only means we choose 2π ≤ Θ i,k < 4π. Now that a ik + N ik and b ik have the same parity, so do x ik and y 0ik . This proves the claim that y 0ik = x ik (mod 2) is a physical requirement for properly chosen Θ i,k and Θ i0,k . With the set {v i } that are linked to unit flux 2π N k e k , the remaining three-loop braiding statistics are easy to reconstruct. To do that, we use the following general properties of Abelian three-loop braiding statistics: The relation (A8a) follows immediately from the definition of exchange statistics. The relation (A8b) comes from the fact that braiding α around β is topologically equivalent to braiding β around α, while both are linked to γ. Equations (A8c)-(A8f) will be referred to as linearity reations. They follow from the fact that braiding and exchanging of loops commute with the two fusion processes of loops, depicted in Fig. 3. More discussions on these linearity relations can be found in Refs. 2 and 7. These works only consider the case that charge excitations are bosonic. Nevertheless, the linearity relations hold regardless of the exchange statistics of charge excitations.
We can now use the two types of fusions in Fig. 3 and the linearity relations (A8c)-(A8f) to obtain braiding statistics between vortices that carry general gauge flux. Also, one can attach charge the vortices to exhaust vortices that carry the same gauge flux. Accordingly, the full set of three-loop braiding statistics indeed can be reconstructed out of the topological invariants. Hence, they are equivalent.

Appendix B: Proofs of the Constraints
In this appendix, we prove the constraints (8a)-(8i) through various physical arguments, except that Eq. (8h) remains a conjecture. The proofs heavily rely on the general properties of Abelian loop braiding statistics Eqs. (A8a)-(A8f).
Proofs of Eqs. (8a), (8b), (8c) and (8g).-First of all, Eq. (8a) follows immediately from the relation (A8a) and the definitions of Θ ii,k and Θ i,k . The constraint Eq. (8b) follows from the relation (A8b). The constraints Eqs. (8c) and (8g) only involve mutual braiding statistics between loops. The fact that there exist fermionic charge excitations does not matter for mutual statistics. Accordingly, they can be established using exactly the same arguments as those given in Ref. 7.
Proof of Eq. (8d).-Consider the thought experiment shown in Fig. 4, where we have N k identical copies of a γ base loop linked with α and α ′ loops. Here, α and α ′ are the same type of loops; we put a prime on the latter just for notational distinction. The loops carry gauge flux φ α = φ α ′ = 2π Ni e i and φ γ = 2π N k e k . Now imagine that we exchange α and α ′ in each copy. The total Berry phase is obviously given by N k θ α,γ . Then, we perform a two-step deformation on the exchange process (Fig. 4): we first fuse the α (α ′ ) loops into a bigger A (A ′ ) loop that is linked with all the γ base loops; then we fuse all the γ loops together to a C loop. It is not hard to see that C carries no gauge flux. Therefore, the original exchange process is deformed to a process of exchanging two loops A and A ′ that are not linked to any base loop (Note that A and A ′ are the same type of loops). Applying the linearity relation (A8f) to the above deformation process, we obtain where the right-hand side is the statistical phase associated with exchanging A and A ′ , and q A is the gauge charge carried by the unlinked loop A.
Proof of Eq. (8e).-Consider a base loop γ that is linked with N i copies of α loops. The loops carry gauge flux φ α = 2π Ni e i and φ γ = 2π N k e k . Using the relations (A8d) and (A8a), it is not hard to see that the exchange statistics of the α loops as a whole is given by N 2 i θ α,γ , which equals N i Θ i,k . On the other hand, fusing the α loops together gives a gauge charge q, whose exchange statistics is πq 0 . Accordingly, we should have In addition, we consider N0 2 copies of β loops that are also linked to γ. The β loops carry gauge flux φ β = 2π N0 e 0 . According to the Aharonov-Bohm law, the mutual statistics between q and the N0 2 copies of β is πq 0 . That is, the mutual braiding statistics between N i copies of α as a whole and N0 2 copies of β as a whole is πq 0 . Using the linearity relations, we find the latter is also given by N0Ni 2 θ αβ,γ . Therefore, we have where the second equality holds when N i is even.
The second thought experiment.
thought experiment in Fig. 4, but now with N ik copies of the linked loops α, α ′ and γ. Similarly to Eq. (B1), we obtain where we assume N ik is even. Note that this q A is different from that in Eq. (B1), since we start with a different number of copies of linked loops in Fig. 4. Next, we consider another thought experiment, shown in Fig. 5. We start with the same N ik copies of linked loops α, α ′ and γ as in Fig. 4. Then, we fuse all the α loops to a big loop A, but we shrink the α ′ loops and fuse them onto the γ loops such that γ turns to γ ′ . The gauge flux carried by γ ′ is the same as that of γ, i.e., φ γ ′ = φ γ . Next, we fuse the N ik copies of γ ′ together and obtain an excitation C ′ . It is not hard to see that C ′ carries no gauge flux, hence it is actually a charge excitation.
To proceed, we create a pair of loops α ′′ andᾱ ′′ , with the loop α ′′ carrying unit flux of 2π Ni e i . We imagine braiding α ′′ around C ′ . Since C ′ is a pure charge, the statistical phase is given by the Aharonov-Bohm phase 2π Ni q C ′ ·e i . On the other hand, since C ′ is composed of N ik copies of loops γ ′ , the linearity relations tell us that the braiding statistical phase is also given by N ik θ α ′′ γ ′ ,A = Θ ik,i . Therefore we obtain We now make use of Eqs. (B4) and (B5) to show the constraint Eq. (8f). Recall that we start with N ik identical α-α ′ -γ links in both Figs. 4 and 5. Each link should carry a well-defined overall charge, which we denote as q link . Because of charge conversation, the overall gauge charge carried by the excitations does not vary in any step of the thought experiments. Accordingly, we have Adding together Eqs. (B4) abd (B5), and using Eq. (B6), we arrive at where we have used the fact that q link is an integer vector.
Finally, we argue that A must be a bosonic charge, i.e., πq A · e 0 = 0. To see that, consider the exchange statistics of C ′ in Fig. 5: Accordingly, C ′ is bosonic. Then, using Eq. (B6) with the assumption that N ik is even, we immediately conclude that A is also bosonic. Hence, we prove the constraint Eq. (8f). Proof of Eq. (8i). -Finally, we argue for the last constraint Eq. (8i). Physically, it follows from the requirement that the flux loops can not have any chiral modes. This is also new to fermionic theories (for bosonic ones, by condensing the bosonic gauge charges on the flux loops one can always make them gapped), and imposes a nontrivial constraint on the braiding statistics of the 2D topological phase obtained from dimensional reduction. We recall the following bulk-boundary relation in a 2D (bosonic) topological phase 57 : Here D = a d 2 a is the total quantum dimension of the 2D theory, a runs over all types of anyons, d a and θ a are the quantum dimension and topological twist of an anyon of type a, respectively. The quantity c − is the chiral central charge, and c − = 0 for nonchiral theories.
Let us apply Eq. (B9). We fix a base loop, and dimensionally reduce the 3D gauge theory. The anyon types in the 2D theory can be labeled by a tuple (q, m), where q labels the gauge charges, and φ i = 2π Ni m i labels the gauge fluxes. These anyon excitations correspond to the vortex loops that are linked to the base loop in the original 3D theory. Assuming all the 2D anyons are Abelian, we have d a = 1 and D = |G| = K i=0 N i . The topological twist of (q, m) is given by where we have used α i to denote the loop that corresponds to the anyon (0, e i ) after dimensional reduction, and have taken the base loop to carry gauge flux 2π N k e k . Next we insert the expression of θ q,m into Eq. (B9) with the requirement c − = 0, and perform the summation over q and m. We first sum over q, and the relevant part is Combining this part with the rest, we have According to the definitions of the topological invariants, the right-hand side of the above equation is given by . (B13) Putting together Eqs. (B12) and Eq. (B13), and properly rewriting the expressions, we obtain the constraint Eq. (8i).
Here B is a bilinear form on A: B(x + y, z) = B(x, z) + B(y, z) and B(z, x + y) = B(z, x) + B(z, y). In our case, we have R a,b = e 2πiB(a,b) , where B(a, b) ∈ Q/Z. Notice that we do not necessarily have B(x, y) = B(y, x). For most of the calculations, we actually have B(y, x) = −B(x, y) mod Z.
The most important property of the cup product is Therefore, if δf = δg = 0, δ(f ∪ g) = 0. One can actually show that the cup product defines a product of cohomology classes. Cup product to some extent is the discrete version of wedge product of differential forms. We also define higher cup product 58 . For our purpose, we mostly just need ∪ 1 : . . , j, j + q, · · · , p + q − 1), g(j, . . . , j + q)].
They satisfy the property: Notice that the sign on the LHS is reversed compared to the usual formula (e.g. see Ref. 58) due to the skew symmetry of the bilinear form B.
Generally, higher cup products satisfy: We list explicit expressions for the ∪ 1 that we will use: Let us define so that We simplify the Pachner move equation for Z + 0 and Z + 1 individually first: In the last step we use the 3-cocycle condition of β.
Combining the two pieces, we get Since β is valued in the transparent center, R β0123,f345 R f345,β0123 = 1 for all f . Therefore the invariance under Pachner moves requires that Appendix E: Gauge Invariance of the Partition Function We will show below that the action can be thought as a topological gauge theory for a 2-form gauge field f and a 1-form gauge field g. We will only consider the case with λ = 0.
To show that the f ijk 's can be regarded as a 2-form gauge field, we need to show that the action is invariant under a 1-form gauge transformation f → f + δξ where ξ is a 1-cochain. The variation of the action is So the partition function does not change.
To establish that g ij 's are 1-form gauge fields requires more work. Under gauge transformations g → g + δh where h is a 0-cochain. First one needs to preserve the flatness condition, so the gauge transformations also affect the 2-form gauge fields. Due to the 3-cocycle condition of β, we can write Here ζ is a 2-cochain. Explicit expressions of ζ can be obtained, but extremely tedious in the general case. So we will illustrate by performing a gauge transformation on a single vertex i: now g ij → g ij + h i while the others remain unchanges. The 3-cocycle β(g ij , g jk , g kl ) transforms as So we can set ζ imn = −β(h i , g im , g jn ).
To preserve the flatness condition δf = β, f has to transform as Illustration of a prism.
Using the formula Eq. (C5), we can write Neglecting the boundary term, Eq. (E4) becomes Now applying the formula Eq. (C6), we obtain up to a boundary term.

Appendix F: Dimensional Reduction
First we need to choose a triangulation of M 3 × S 1 . We start from an open manifold M 3 × D 1 (where D 1 stands for an interval), and triangulate the two boundaries into 3-simplicies denoted by [ijkl] and [i ′ j ′ k ′ l ′ ], respectively. We will exploit an ordering in which i < i ′ for all vertices i in the simplicial triangulation of M 3 . We then identify the two boundaries, i.e. i is identified with i ′ , to obtain M 3 × S 1 . The G flux through S 1 is measured by a Wilson loop along S 1 , which is the G label on ii ′ (it is easy to see that all jj ′ should have the same label by the flatness condition). We will set g ii ′ = h from now on.
The basic building block of the triangulation is a prism Next, we divide the 2-simplices into two types: the "in-plane" ones, which lie entirely inside the 3-manifold M 3 , i.e. f ijk ≡ f i ′ j ′ k ′ , and those going "out of plane", e.g. f ii ′ j ′ . They need to be treated differently. This is already evident when we examine the twisted flatness conditions: the flatness conditions for the in-plane fields are essentially properties of M 3 , and may depend on the topology of the manifold. The flatness conditions for the out-of-plane fields can be dealt with explicitly, which we will analyze now.
Let us consider the twisted flatness conditions on a 3D prism [ijki ′ j ′ k ′ ], which is further triangulated into three 3-simplices [ijkk ′ ], [ijj ′ k ′ ] and [ii ′ j ′ k ′ ]. Let us write the flatness conditions out: The meaning of these conditions is uncovered by considering (F1a)−(F1a)+(F1c): Here m ij are defined as m ij = f ii ′ j ′ − f ijj ′ , and n(g ij , g jk ) is given by (F3) Here i h β is called the slant product of β: Now we need to evaluate the partition function. First we collect the contributions from T's on a single prism We can use flatness conditions on [00 ′ 2 ′ 3 ′ ] and [0133 ′ ]: to rewrite Let us collect the following factors involving f ijk 's: To see that these terms do not contribute to the partition function, we recall the following property of the cup product: Recall that δf = β, δm = −n, Eq. (F8) implies Here ∼ means up to a boundary term. To further simplify the expressions let us do the following gauge-fixing: we fix all f ijj ′ = 0, by using the gauge degrees of freedom on ij ′ , and then m ij = f ii ′ j ′ . So f 00 ′ 3 ′ − f 00 ′ 2 ′ = m 03 − m 02 = m 23 + n 023 , f 033 ′ − f 133 ′ = 0. Notice that after the gauge-fixing m ij , we also need to multiply the partition function by a factor of |A| to correctly normalize it. In total the partition function should be multiplied by |A| |∆1(M3)| . Now we have obtained the following expression for the partition function on a prism: i h ω(g 01 , g 12 , g 23 ). (F10) We separate the weights associated to m and f : where α(g 01 , g 12 , g 23 ) = R n012+n023,β0123 R n023−β(h,g02,g23),n012 R β(g01,g13,h),n123 i h ω(g 01 , g 12 , g 23 ). (F12) The exact expression of α is not important for our analysis.
Let us now take care of the normalization factors. We have Thus the normalization factor becomes To summarize, we have found that the partition function on M 3 × S 1 can be written as It should be clear that in this expression ∆ k ≡ ∆ k (M 3 ). In the following we discuss the physical interpretation of the sum over the "in-plane" fields f . We will only consider the simple case ξ = 0. The sum over f can be evaluated: In the following we will write Let us calculate the number of 1-and 2-cocycles: where R a,b = (−1) ab , β is a 3-cocycle in H 3 [G, T ], and λ is a 2-cocycle in H 2 [G, T ]. When ν is not a 5-coboundary, we say that the corresponding model has an H 5 [G, U(1)] obstruction. The main purpose of this appendix is to determine when the twisted Crane-Yetter models are obstruction-free, i.e. when the ν in Eq. (G1) is a 5-coboundary. We only discuss the case that G is an Abelian group and T = Z N0 .
To do that, we define six quantities {Θ i,l,m , Θ ij,l,m , Θ ijk,l,m , Ω i , Ω ik , Ω ijk } for a general 5-cocycle ν ∈ H 5 (G, U(1)), where G is any finite Abelian group Z Ni . These quantities have the important property that they are defined in a way such that they are invariant under a coboundary transformation ν → νδµ. Hence, we call these quantities invariants for H 5 [G, U(1)]. We claim that a 5-cocycle is a 5-coboundary if and only if all the six corresponding invariants vanish. We define these invariants for general 5-cocycles in appendix G 1, and we prove the claim in appendix G 2. Finally, we apply the invariants to the specific 5-cocycle given in Eq. (G1).

Defining the invariants
To define the invariants, let us first define the following functions: i a ν(g, h, k, l) = ν(a, g, h, k, l)ν(g, h, a, k, l)ν(g, h, k, l, a) ν(g, a, h, k, l)ν(g, h, k, a, l) (G2) The function i a ν is usually called the "slant product" of ν. It is actually a 4-cocycle in H 4 [G, U(1)], when a is treated a parameter. One may continue to apply the slant product on i a ν: i b i a ν(g, h, k) = i a ν(g, b, h, k)i a ν(g, h, k, b) i a ν(b, g, h, k)i a ν(g, h, b, k) where i b i a ν is a 3-cocycle and i c i b i a ν is a 2-cocycle. In addition, we also define the following function i a,b ν(g, h, k) = ν(g, h, k, a, b)ν(g, a, h, k, b) ν(g, h, a, k, b)ν(a, g, h, k, b) ν(g, h, a, b, k)ν(a, g, h, b, k) ν(g, a, h, b, k) ν(g, a, b, h, k) ν(a, g, b, h, k) ν(a, b, g, h, k) (G4) The function i a,b ν, however, is not a 3-cocycle.
With these functions, we now define the invariants for H 5 [G, U(1)] for Abelian group G = i Z Ni . Let e i be the generator associated with the Z Ni subgroup of G. First, we define the following invariants for a given 5-cocycle ν: where the products are taken over all prime numbers p. Then, we have the following group isomorphisms Z Ni = Z 2 r 2 × Z 3 r 3 × Z 5 r 5 × · · · Z Nj = Z 2 s 2 × Z 3 s 3 × Z 5 s 5 × · · · Z N k = Z 2 t 2 × Z 3 t 3 × Z 5 t 5 × · · · Let e p i ≡ Ni p rp e i be the generator associated with the Z p rp subgroup in Z Ni , and e p j , e p k are similarly defined. In the case that r p ≤ s p ≤ t p , we define i e p k (e p i + e p j , ne p i + ne p j , e p i + e p j , me p i + me p j ) i e p k (e p i , ne p i , e p i , me p i )i e p k (e p j , ne p j , e p j , me p j ) If r p , s p , t p are in different orders, Ω p ijk are defined similarly with a corresponding permutation of indices i, j, k in Eq. (G10). At the end, we define the total invariant Again, one can show that Ω ijk is invariant under a coboundary transformation ν → νδµ.

Completeness of the invariants
Let us now show that the invariants {Θ i,l,m , Θ ij,l,m , Θ ijk,l,m , Ω i , Ω ik , Ω ijk } are complete, in the sense that they have the resolution to distinguish every cohomology class in H 5 [G, U(1)]. To do that, we perform a counting argument. First of all, for Abelian group G = i Z Ni , the cohomology group is given by That means the invariants can take at most |H 5 [G, U(1)]| distinct values. If we are able to show that the invariants can take exactly |H 5 [G, U(1)]| distinct values, we prove the invariants are complete.
To do that, we evaluate the values of the invariants for the following explicit 5-cocycles ν 1 (a, b, c, d, e) = exp where P ijk , Q ijkl , R ijklp are integer parameters. For simplicity, we assume R ijklm = 0 if any of its two indices are equal. We have used integer vectors a = (a 1 , a 2 , . . . ) to denote the group elements of G with 0 ≤ a i < N i , and [b j + c j ] is defined as b j + c j modulo N j . We use additive convention for group multiplication of the Abelian group G. One may check that for ν 1 and ν 2 are indeed 5-cocycles. Inserting the expression of ν 2 into the definition of Θ ijk,l,m , Θ ij,l,m , Θ i,l,m , we find that Θ ijk,l,m = − 2π N ijklm σ sgn (σ)R σ(i)σ(j)σ(k)σ(l)σ(m) (Note that Ω i , Ω ij , Ω ijk can be evaluated but the expressions are complicated, so we do not list them here. This does not affect the counting argument below. ) Insertiging the expression of ν 1 in to the definitions of Ω i , Ω ik , Ω ijk , we find that Θ ijk,l,m = Θ ij,l,m = Θ i,l,m = 0 and Ω i = 2π P iii N i Ω ik = 2π P iik + P iki + P kii N ik Ω ijk = 2π P ijk + P ikj + P kij + P jik + P jki + P kji N ijk X (G15) where X is an integer with the property that X and N ijk are coprime.
Let us count how many distinct values these invariants can take. First, Θ ijk,l,p is fully antisymmetric, and it can take N ijklp different values. One can show that N ijlp Θ ij,l,p = 0, and it is symmetric in i, j and antisymmetric in l, p. A more careful calculation shows that for fixed indices i = j = l = p, Θ ij,l,p and those related to Θ ij,l,p by index permutations can take N 3 ijlp distinct values. The invariant Θ i,l,p is antisymmetric in l, p. For fixed indices i = j = p, Θ i,l,p and those related by index permutations can take N 3 ilp distinct values. Hence, Θ ijk,l,p , Θ ij,l,p , Θ i,l,p together can take N Θ different values with Next, we count the possible number of values for Ω i , Ω ij , Ω ijk from ν 1 . Since for ν 1 we have Θ ijk,l,m = Θ ij,l,m = Θ i,l,m = 0, these possible values of invariants are distinct from those from ν 2 . The invariant Ω i obviously can take N i distinct values. For given i = k, Ω ik and Ω ki are independent. They together can take N 2 ik distinct values. The invariant Ω ijk is fully symmetric and it can take N ijk distinct values (note that in the definition of Ω ijk , it is only symmetric in i and j. This full symmetry in all three indices is only a consequence of the specific cocycle ν 1 ). Accordingly, we have that the invariants Ω i , Ω ij , Ω ijk can take N Ω distinct values with Putting all together, the invariants can take N Θ N Ω = |H 5 [G, U(1)]| distinct values in total. Hence, they are complete.

Evaluating obstruction
We now evaluate the values of the invariants {Θ i,l,m , Θ ij,l,m , Θ ijk,l,m , Ω i , Ω ik , Ω ijk } for the 5-cocycle given in Eq. (G1). That 5-cocycle depends on a 3-cocycle β in H 3 [G, Z N0 ] and a 2-cocycle λ in H 2 [G, Z N0 ]. Below, we work for the following explicit β and λ: where p ij , p ijk , q i are integer parameters, and for simplicity, we assume that p ijk = 0 if any two of the indices are the same.
We now insert the above expressions of β and λ into the expression of ν in Eq. (G1), and further insert ν into the definitions of the invariants {Θ i,l,m , Θ ij,l,m , Θ ijk,l,m , Ω i , Ω ik , Ω ijk }. After a long tedious calculation, we find that (1 + q i + q j )p ijk + π N 0 N ijk0 (q ipjk + q jpki + q kpij ), (when i = j = k) Θ ijk,l,m = 0 (G19) wherep ijk ≡ p ijk + p jki + p kij − p kji − p jik − p ikj andp ij = p ij + p ji . (Note that we only calculated Ω ik and Ω ijk with i = j = k for simplicity. For Ω ijk , k is the index such that N k contains the most the factor 2 among N i , N j , N k .) The twisted Crane-Yetter models are obstruction free if and only if all the above expressions vanish.