Dynamics of a lattice 2-group gauge theory model

We study a simple lattice model with local symmetry, whose construction is based on a crossed module of finite groups. Its dynamical degrees of freedom are associated both to links and faces of a four-dimensional lattice. In special limits the discussed model reduces to certain known topological quantum field theories. In this work we focus on its dynamics, which we study both analytically and using Monte Carlo simulations. We prove a factorization theorem which reduces computation of correlation functions of local observables to known, simpler models. This, combined with standard Krammers-Wannier type dualities, allows us to propose a detailed phase diagram, which form is then confirmed in numerical simulations. We describe also topological charges present in the model, its symmetries and symmetry breaking patterns. The corresponding order parameters are the Polyakov loop and its generalization, which we call a Polyakov surface. The latter is particularly interesting, as it is beyond the scope of the factorization theorem. As shown by the numerical results, expectation value of Polyakov surface may serve to detects all phase transitions and is sensitive to a value of the topological charge.


Introduction
Higher gauge theories are physical models which generalize conventional gauge theory by associating degrees of freedom to geometric objects of dimension higher than one. Perhaps the best known example is the p-form electrodynamics [1], whose discretized version, a natural generalization of the Wilson formulation of ordinary gauge theory [2,3], can be expressed in terms of degrees of freedom associated to p-cells, e.g. plaquettes for p = 2. These degrees of freedom are subject to redundancy described by group valued functions on the set of (p − 1)cells. In the case of p = 1 this reduces to degrees of freedom on links with gauge transformations given by arbitrary functions defined on lattice sites.
Already in [1] it was argued that gauge theories with p ≥ 2 are necessarily abelian, essentially because there exist no well-behaved orderings on surfaces. There exists a way to bypass this argument, inspired by higher category theory [4,5,6]. For p not exceeding 2, it is typically formulated in terms of 2-groups [7] or, equivalently, crossed modules [8]. Surface observables in 2-group gauge theories are still valued in an abelian group, but in general they are computed in terms of genuinely non-abelian local degrees of freedom associated to links and plaquettes.
There exists also a concept of (global) higher form symmetries [9], whose relation with higher gauge theories is similar to the relation between ordinary symmetries and gauge theories. Examples of models admitting higher symmetries have been known for a long time, and among gauge theories they are in fact the rule rather than an exception. Nevertheless, systematic study of higher symmetries seems to have begun relatively recently.
Higher gauge theories have been proposed [10,11] as effective field theories describing vacua of conventional gauge theories. They also provide interesting examples [12,13,14,15,16,17] of Topological Quantum Field Theories (TQFTs) [18,19], and hence are expected to describe certain gapped topological phases of many body quantum systems. In [12] the existence of Symmetry Protected Topological (SPT) phases protected by higher symmetries was proposed. Another motivation to study higher gauge theories is provided by its relation with bosonization in arbitrary dimension [20,21,22]. Furthermore, certain models in string theory may be described as higher gauge theories [23].
Yetter's model [13] is a TQFT based on a crossed module of finite groups. Its hamiltonian realizations resembling the Kitaev's toric code were constructed in [24,25,26]. In [27] a common generalization of the Yetter's model, 2-form Z n electrodynamics and lattice Yang-Mills theory has been proposed. It is a genuinely dynamical model, formulated in the hamiltonian formalism, which reduces to a TQFT only in certain limits. In this work we consider an analogous model formulated in terms of state sums (or discrete functional integrals). We focus on one relatively simple crossed module, but some of our results are true in general. In order to make the paper more accessible, we have decided to define everything explicitly using notations standard in lattice gauge theory. We refer to [27] for an exposition of the slightly more involved formalism of crossed modules and proofs of various algebraic facts used in this text.
Full definition of the model under consideration is given in subsection 2.1. Its extended observables are discussed in subsection 2.2. We identify topological charges and higher symmetries: 1-form symmetry Z (1) 2 and 2-form symmetry Z (2) 2 . We discuss the theoretical possibility of symmetry breaking and provide suitable order parameters. In subsection 2.3 we show that computation of a large class of observables, including all local observables, may be reduced to calculation of averages in simpler models: 1-form Z 2 gauge theory and 2-form Z 2 gauge theory. This includes the statement that plaquette observables (constructed from link degrees of freedom in the usual way) are uncorrelated with cube observables (constructed from plaquette degrees of freedom), which is not obvious from the form of the action. This factorization theorem is not valid for the surface observable which is the order parameter of the Z (2) 2 symmetry. In subsection 2.4 we use the factorization theorem and Kramers-Wannier type dualities to formulate a proposal for the phase diagram. We describe critical points, symmetry breaking patterns and renormalization group fixed points governing the infrared physics.
Section 3 is devoted to Monte Carlo study of the proposed model in dimension D = 4. Simulation algorithm is described in subsection 3.1. Since the general model is fairly standard, we discuss in detail only those aspects that are specific to the case at hand. In subsection 3.2 we present numerical results for expectation values of local observables. These results confirm Figure 1: Independent degrees of freedom are associated to links m µ (x) and faces n µν (x). The latter should not be confused with plaquette observables f µν (x) constructed from links.
the phase structure obtained from duality arguments. The most interesting, in our view, results of simulations are presented in subsection 3.3. They concern behaviour of order parameters for higher symmetries Z 2 and Z 2 . It is found that order parameters for the latter not only exhibit sharp dependence on both coupling constants, but also are sensitive to the topological charge.
Finally, in the appendix A, we discuss construction of non-spherical surface observables using the general language of crossed modules.
2 Description of the model

Degrees of freedom, action and gauge freedom
Coordinates of a lattice site form a tuple x = (x 1 , ..., x D ) with integer x µ . Unit vector in the direction µ ∈ {1, ..., D} will be denoted by µ. We choose periodic boundary conditions, i.e. x µ is identified with x µ + L µ , where L µ is the spatial extent of the system in the µ-th direction.
Addition and multiplication in Z 4 = {0, 1, 2, 3} is always performed modulo four. We consider a model with degrees of freedom of two types, both valued in Z 4 (see Fig. 1): • m µ (x), associated with the link between x and x + µ, • n µν (x), associated with the square with corners x, x + µ, x + µ + ν, x + ν (called face).
They are subject to a constraint (for every x and µ < ν) called fake flatness: (2.1) The right hand side is a plaquette built of m variables as in the ordinary lattice gauge theory. It is convenient to denote it by f µν (x). We note that n µν (x) is determined by the link variables only modulo two and that f µν (x) has to be even (but m µ (x) not necessarily so).
Out of elementary degrees of freedom one may construct observables associated to cubes: The six terms in this formula correspond to six sides of a cube. It can be shown that fake flatness enforces all g µνρ (x) to be even.
Observable f µν (x) is the Wilson line along the boundary of an elementary rectangle. In the present model it is possible to construct also higher dimensional analogues of Wilson lines, which could be called Wilson surfaces. Observable g µνρ (x) is the Wilson surface along an elementary cube.
We choose the following action: Degrees of freedom superficially seem to interact with each other, since they are related by the fake flatness condition and since g µνρ (x) (which enters the action directly) depends on both degrees of freedom. However, as it will be demonstrated later, this interaction does not affect local dynamics, i.e. plaquettes are uncorrelated with cubes and furthermore correlation functions of plaquettes and cubes depend only on J 1 and only on J 2 , respectively. On the other hand, the impact of the interaction can be seen in averages of nonlocal order parameters. Numerical evidence supporting this statement is presented in section 3.
The fake flatness constraint (2.1) and the action (2.3) are invariant under gauge transformations of two types. Gauge transformations associated to sites are parametrized by elements ξ(x) ∈ Z 4 . They act according to the formulas: Gauge transformations associated to links are parametrized by ψ µ (x) ∈ {0, 2} Z 4 and act as (2.5b) Only gauge invariant quantities will be regarded as observables. In this work we consider f µν (x), g µνρ (x) and order parameters described in subsection 2.2.

Nonlocal order parameters and symmetries
Polyakov loop, a particular Wilson line winding around one of the directions of the lattice, is defined by the formula Its possible values are ±1 and ±i, in contrast to plaquette observables which take only two possible values. For configurations with f µν (x) = 0 for all x, value of p µ (x) is independent of x. This is not true for general configurations. On the other hand, quantity Q µ defined by is independent of x, which follows from the fact that all f µν (x) are even. We will call it the topological charge. Each Q µ may take two possible values, 1 or −1, so the whole set of field configurations decomposes into 2 D disjoint sectors. We note that local constraint-preserving transformations in the set of all field configurations cannot change the topological charge, since that requires changing p µ (x) for all x.
For every µ there exists a symmetry of the action which leaves all plaquettes, cubes and {p ν (x)} ν =µ unchanged, but flips the sign of p µ (x) (and hence preserves Q µ ). It is given by (2.8b) We will call it the electric 1-form symmetry. As a consequence of this symmetry the expectation value of p µ (x) vanishes.
We are unaware of a symmetry which changes p µ (x) by a factor of i (and hence flips the sign of Q µ ). Nevertheless, it will turn out to be useful to consider the Q µ -reversing transformation (2.9b) It preserves fake flatness and all f µν (x), so it is a symmetry of S 1 . However, it changes values of cube observables, so it is not a symmetry of the full action.
We would like to address the question whether the symmetry (2.8) can be spontaneously broken. We insist on gauge invariance and locality of the action, so it is not possible to include a symmetry breaking term in the action. On the other hand, in a putative phase with unbroken symmetry the infinite volume limit of the volume average of p µ (x) in a fixed typical field configuration is expected to vanish. This happens for small J 1 , because then plaquette observables fluctuate strongly, so the signs of p µ (x) and p µ (y) are essentially independent if the transverse distance |x − y| ⊥ = ν =µ (x ν − y ν ) 2 is large. More precisely, p µ (x)p µ (y) −1 may be understood as a Wilson loop bounding area L µ |x − y| ⊥ , so its expectation value is expected to decay exponentially with |x − y| ⊥ .
To quantify the above discussion, we consider with the sum taken over x in a plane transverse to the µ-th direction. Here V ⊥ = ν =µ L ν is the transverse volume. Squaring this definition we find There are V 2 ⊥ terms in this sum, each of which has modulus one. After taking expectation value, only O(V ⊥ ) terms, with |x − y| ⊥ comparable to the correlation length survive. Therefore the By spontaneous breaking of the symmetry (2.8) we shall understand violation of this scaling law, so that P µ remains nonzero in the limit of infinite transverse volume: Note that it may still be true that P µ → 0 as L µ → ∞.
There exists a surface observable analogous to the Polyakov loop. It may be thought of as a Wilson surface winding around two lattice directions. Its construction is slightly more involved. We choose a plane through a fixed site x parallel to directions µ < ν. Morally speaking, we would like to add n µν (y) with y running through all sites in the chosen plane. However, this does not give a gauge invariant quantity. To fix this, we have to choose for every y a path from y to x (which we take to be entirely contained in the chosen plane) and weigh n µν (y) by a parallel transport factor (−1) mρ(z) , where the product is taken over all links forming the chosen path. Then the sum, denoted Σ µν (x), is gauge-invariant and even. This is discussed in the broader context of crossed modules in the appendix A. We define 14) which will be called the Polyakov plane.
We remark that our notation is fully justified only if either Q µ = Q ν = 1 or all f µν (x) vanish, because otherwise p µν (x) depends on the choice of paths needed to construct it. This hints at the possibility that expectation values of p µν (x) may depend both on the two coupling constants and on the topological charge. This will be corroborated by results in section 3.
There exists a symmetry which flips the sign of p µν (x): We will call it electric 2-form symmetry. It implies that the expectation value of p µν (x) vanishes.
By analogy with the Polyakov loop, we consider the quantity where V ⊥ = ρ =µ,ν L ρ and the sum is taken over a plane transverse to µ and ν. We will say that the symmetry (2.15) is broken if P µν has nonzero limit as V ⊥ → ∞.

Reduction of dynamics to simpler models
In this subsection we will show how to express certain averages with respect to the action (2.3) in terms of averages in simpler models. We will make use of constraint-preserving moves in the space of field configurations. Firstly, the link moves: with arbitrary ψ µ (x) ∈ Z 4 . They reduce to gauge transformations (2.5) if ψ µ (x) is even. In general they change the value of plaquette observables f µν (x), but not of the cube observables g µνρ (x). Secondly, the face moves: . They preserve f µν (x), but change values of g µνρ (x).
Moves described above generate a group. Every move may be represented as a sequence of local moves with only one nonzero ψ µ (x) or χ µν (x). Since m µ (x) are always either unchanged or shifted by an even amount, topological charges Q µ are invariant.
We claim that any two configurations with equal topological charges can be related by a sequence of local moves and a gauge transformation. Indeed, first consider two configurations with equal m µ (x). Then, by fake flatness, all n µν (x) differ by even numbers, so the two configurations are related by a face transformation. This reduces the proof of the claim to showing that m µ (x) can be made equal by a sequence of link moves and a gauge transformation. The only gauge invariant functions of m µ (x) are Wilson lines, which can be taken along contractible loops or non-contractible loops. The former are expressible in terms of f µν (x) and have to be even. The latter are also even on the account of the assumption about topological charges, since every loop can be built of contractible loops and Polyakov loops. This proves that up to pure gauge terms, the difference of m µ (x) is even. Thus they are related by a transformation of the form (2.17a).

The average of an observable
in which the sum over m, n is restricted by the constraint. Function f (Q) is a weight given to the sector with topological charge Q. The simplest choice is f (Q) = 1, while restriction to Q = Q with fixed Q corresponds to f (Q) = δ Q,Q . We consider an observable of the form • O 1 can by expressed solely in terms of plaquette observables f µν (x) (thus it can be P µ ), • O 2 is invariant with respect to gauge transformations and link moves, e.g. it is an arbitrary function of cube observables g µνρ (x).
We define the quantity It is invariant with respect to gauge transformations and link moves of m variables, so it depends on m only through Q µ . Therefore we write We divide the summation over m into topological sectors. The sum over m with fixed Q will be denoted by index m|Q: does not depend on Q by symmetry (2.9) of S 1 . Finally: In the remaining sum over m we have configurations of m variables such that every Wilson loops is even. Such configuration is gauge equivalent to one with all m µ (x) even. Furthermore, every gauge orbit has 4 N 1 −1 representatives (in which N 1 is the number of links), out of which 2 N 1 −1 is such that all m µ (x) are even. Therefore we may restrict the sum over m to configurations with even m µ (x) at the small cost of including a factor 2 N 1 −1 . Then the sum over m gives the Wegner model [28], so This gives from which we can draw the following conclusions: • O 1 does not depend on J 2 and weights f , and is equal to the average in Wegner's model, This factorization theorem is the main result of this section. We would like to remark that its derivation remains valid also for models based on general crossed modules of finite groups, as long as the action is a sum of a term depending only on plaquette observables and a term depending only on cube observables. This observation follows from the fact that the presented proof relies only on general properties of gauge transformations and constraint-preserving moves. These were discussed in detail in [27].
Next we turn to the question on how O 2 depends on the topological charge sector. We will argue that for thermodynamic quantities dependence becomes negligible in the infinite volume limit. This will be confirmed already for quite small lattices by results of simulations presented in section 3.
Consider, for concreteness, the case Q 1 = −1, Q µ = 1 for µ = 1. Such choice of topological charge may be realized by the gauge field m µ (x) = δ µ,1 δ xµ,0 . (2.26) It is supported on a plane, so switching it on (without modifying n µν (x) variables) may change the value of at most D−1 L µ cubes. Hence we have It follows that obeys an estimate Taking logarithms gives an estimate on the difference of free energies per unit volume: We recall that the free energy is an extensive quantity. On the other hand, the right hand side divided by the volume decays as L −1 1 as L 1 → ∞. We conclude that in the infinite volume limit, the free energies per unit volume become equal in all topological sectors.

Phase diagram proposal for D = 4
In this subsection we restrict attention to dimension D = 4, although some parts of the discussion are valid also for other dimensions.
In the case D = 4, Wegner's model has a single phase transition [28], which is of first order. Its exact position may be calculated using Kramers-Wannier type self-duality 1 [29]. There exist two renormalization group fixed points at J 1 = 0 and J 1 = ∞. Two phases may be interpreted as their basins of attraction.
At the point J 1 = 0, degrees of freedom become completely random and hence the theory is trivial. Effect of a small, but nonzero J 1 may be calculated using the strong coupling expansion [30]. One finds that Wilson loops obey the area law, and hence the electric 1-form symmetry is unbroken.
At J 1 = ∞ the system is constrained to configurations which minimize the action. Thus all plaquette observables vanish and Polyakov loops become independent of position. Up to gauge transformations, minima of the action are labeled by values of Polyakov loops. In Wegner's model each P µ takes 2 possible values, so there exist 16 minima. They all have the same value of the action, because they are connected by the electric 1-form symmetry. However, in order for the system to get from one minimum to another using local moves only, it has to overcome an infinite action barrier. Even for finite J 1 (but large, so that a typical configuration is close to a minimum) the height of the barrier is of order J 1 V ⊥ , so one may expect the electric 1-form symmetry to be broken.
The link variable sector of our model is slightly different in that the Polyakov loop takes four, rather than two possible values. However, it becomes essentially equivalent to the Wegner's model after restricting to a single topological charge sector.
Next we turn to the local dynamics of plaquette degrees of freedom. There exists a Kramers-Wannier duality between W 1,trivial (J 2 ) and the Ising model partition function 2 with sinh(2J 2 ) sinh(2J Ising ) = 1. (2.31) In the Ising model one expects a single continuous phase transition 3 whose position reported in [32] is J crit Ising = 0.149647 (5). This corresponds to a continuous phase transition in our model at The critical point of the Ising model is expected to be described by a massless scalar field theory. It admits one relevant perturbation, given by the mass term. Therefore the fixed point at J 2 = J crit 2 is repulsive. The only other fixed points at J 2 = 0 and J 2 = ∞ describe physics in phases J 2 < J crit 2 and J 2 > J crit 2 , respectively.
Quite analogously to the Wegner's model, the electric 2-form symmetry is unbroken in the small J 2 phase. The situation is much more interesting for large J 2 . To gain some orientation about this case, we consider the limit J 2 = ∞, in which configurations are constrained to minimize S 2 . As shown in the appendix A, Polyakov surfaces p µν (x) become independent of x if in addition either J 1 = ∞ (i.e. for configurations minimizing also S 1 ) or if topological charges are trivial. Therefore we expect that the 2-form symmetry is broken if J 2 > J crit 2 and J 1 > J crit 1 . In the phase J 2 > J crit 2 , J 1 < J crit 1 we can make this conclusion only for the sector with trivial topological charge. On the other hand, numerical results in section 3 show that in the sector with Q µ = −1 the symmetry is restored. We find this result quite striking.
We remark that the four renormalization group fixed points described here may be identified with four integrable hamiltonians described in [27].

Simulation method
In the numerical setup, we keep the extent of three directions equal L 0 = L 1 = L 2 = L, whereas L 3 will be varied separately. We denote the entire volume by V = L 3 × L 3 . We will also use the notation (x 0 , x 1 , x 2 , x 3 ) = (x, y, z, t). (3.1) For any observable we define its statistical expectation value, denoted by · , as the arithmetic mean over samples from a single Markov chain, and in some cases a weighted average of expectation values from multiple Markov chains. In most cases we perform a single simulation where we gather around 10 5 measurements, from which we estimate the average and its standard deviation, taking into account autocorrelations. We do the latter by explicitly calculating the autocorrelation function and integrating it up to the first non-positive element to quantify the autocorrelation time τ int . In the following figures, all data points are shown together with their statistical uncertanties, which however may be smaller than the symbol size and hence invisible. In some cases we have performed up to four parallel simulations in order to increase the statistics and to check for ergodicity.
All simulations are performed using an intertwined application of Metropolis [33,34,35] and over-relaxation steps [36,37,38,39,40]. These are two independent update steps coming in pairs: one for updating the link variables and another to update face variables. We now describe both in more details.
Metropolis steps are based on local changes separately for both kind of degrees of freedom. We use (2.17a) to update the link variables and (2.18a) for the face variables. We remind that by construction such moves preserve the fake-flatness constraint. As a consequence, the move (2.17a) changes both link and face variables. If the constraint was satisfied by the initial configuration, it will be satisfied during the successive application of any of the above changes. Any two configurations can be linked by a finite-length chain of such local movements, which ensures that the simulations are ergodic. Each new configuration ν is obtained from a previous configuration µ by a local change of a randomly chosen degree of freedom and is subject to an accept/reject step with a probability given by p A (µ → ν) = min 1, e S(µ)−S(ν) . (3.2) Over-relaxation steps are made of non-local moves (2.8) and (2.15), which flip the signs of Polyakov lines P µ and Polyakov planes P µν , respectively. Since such transformations do not change the value of the action of a given configuration, they would be always accepted. Hence they are not subject to the accept/reject step. It is known that the incorporation of such moves between Metropolis moves reduces autocorrelation times significantly [36,37,38,39,40].
The local moves (2.17a) and (2.18a) cannot change the value of the topological charge. Hence, the simulation is limited to the topological sector given by the value of the topological charge of the initial configuration. In the following we discuss two independent chains of simulations, one performed in the trivial topological sector (Q µ = 1 for all µ) and the second performed in the sector with Q 0 = −1, see (2.7). We construct the latter by starting from an initial configuration where all the link and face variables are set to 0. Subsequently we set m 0 (0, y, z, t) = 1 for all y, z, t, thus enforcing P 0 = i and hence Q 0 = −1.
The above algorithm with the accept/reject as in (3.2) satisfies the detailed balance condition, which together with the ergodicity of the local moves, guarantees the correctness of the entire algorithm in a given topological sector.
In order to identify the thermalization region of the Markov chain we usually perform an additional simulation with the same parameters, which we start from a so-called hot initial configuration. The latter is constructed by randomizing as much as possible all the degrees of freedom. To be more precise, we set m µ (x, y, z, t) to 0 or 2 with equal probabilities, and subsequently adjust n µν (x, y, z, t) variables to satisfy the fake-flatness constraint. We do this by evaluating all plaquette variables f µν (x, y, z, t) and then setting n µν = 1 2 f µν (x, y, z, t) + q, where q is a random variable taking values 0 and 2 with equal probability. In both simulations, started from a cold and hot configuration, we monitor two local variables (3.3) and (3.4). Recording of relevant observables is started only when the two monitored quantities attain compatible values.

Numerical results for local observables
In this section we discuss two local observables: plaquettes and cubes F = 1 6V x,y,z,t µ<ν f µν (x, y, z, t) , x,y,z,t µ<ν<ρ g µνρ (x, y, z, t) . (3.4) According to the factorization theorem (2.25), we expect that F does not depend on J 2 and G does not depend on J 1 . Our first numerical results confirm these conclusions. In Figure 2 we show the average values F and G as functions of J 1 and J 2 separately. Plots of data obtained in different topological sectors are also indistinguishable, up to statistical uncertainties.
In Figure 2 we demonstrate the dependence of F and G on J 1 (two panels in the upper row) and J 2 (two panels in the lower row) coupling constants. Motivated by our expectations regarding the phase diagram of the system, i.e. existence of four distinct phases, as described in section 2.4, linked to the corners of the phase space given by (J 1 , J 2 ) = (0, 0), (0, ∞), (∞, 0) and (∞, ∞), we select values of J 1 and J 2 representing each phase: We clearly see in Figure 2 that G does not depend on J 1 , i.e. the values are constant and compatible within their statistical uncertainties for the entire range of J 1 values investigated.  Fig. 3 demonstrates that indeed a second order phase transition happens around the expected J crit 2 . In all the cases, results from both topological sectors are shown: Q 0 = 1 and Q 0 = −1, with Q µ = 1 for µ = 0. The data points for the latter are shifted by 0.0025 along the x-axis in order to increase the plot readability.
Similarly, F does not depend on J 2 . Near the location of the expected first order phase transition in J 1 , value of F drops significantly. We demonstrate the nature of this phase transition in the left panel of Figure 3. The panel reproduces the results from [41], where the hysteresis in the average plaquette action in the four-dimensional Wegner model [28] was interpreted as a clear sign of a first order phase transition.
As far as G is concerned, the lower right panel shows a rather smooth dependence. The part of the action proportional to J 2 is a function of G , hence we conclude that also the action itself has a continuous dependence on J 2 . This is in agreement with the expected nature of the phase transition in J 2 being second order. We corroborate this with the results shown in the right panel of Figure 3, where fluctuations of G are shown to exhibit a drastic change around J crit 2 . To be precise, we plot (G − G ) 2 V −1 . Again, all results for G show no dependence on the change in the J 1 coupling constant.
The results shown in Figure 2 provide an illustration of the factorization theorem. Moreover, they support the expected existence of four phases at the four corners of the phase diagram. The more detailed results shown in Figure 3 suggest that the location of the critical couplings where the phase transitions occur, agree within the accuracy of our simulations with the predictions (2.30) and (2.32). Hence, already the simple, local observables such as F and G provide valuable information about the system. We now turn our attention to non-local observables: Polyakov line and Polyakov planes. The latter, as opposed to the former, are not subject to the factorization theorem and hence are expected to have a non-trivial dependence on both J 1 and J 2 .

Numerical results for non-local observables
In this section we discuss results for extended observables. We study in details two such observables: the (volume averaged) Polyakov line, P 0 , winding around the x-direction (2.10) and the Polyakov plane P 01 , winding around the x and y directions (2.16).
Our expectations for these observables in the four possible phases are based on considerations of the system of infinite size in directions perpendicular to the winding directions. We mimic that limit by taking L 3 → ∞, which is the direction perpendicular to both P 0 and P 01 . We discuss our numerical findings below.
We show the numerical results for P 0 and P 01 in Figures 4 and 5 at four pairs of coupling constants as a function of the extent of the lattice in the L 3 direction. In the left panels we gather results obtained at J 1 = 0.43 and J 1 = 0.46 at small J 2 = 0.10, whereas in the right panels we keep the same two values of J 1 but we change J 2 to a large value, J 2 = 1.10. As opposed to the previous section, where F and G were discussed as functions of J 1 and J 2 varying around their critical values, here we study the dependence on the L 3 extent at the four values of coupling constants selected in (3.5) and (3.6).  We start with the discussion of Polyakov lines. The observable P 0 is expected to satisfy the factorization theorem. Indeed, we find that its average value does not depend on the value of the J 2 coupling constant. As a consequence, the left and right panels of Figure 4, showing the results for J 2 = 0.10 and J 2 = 1.10 respectively, look very similar. Two scenarios can be realized as the volume of the lattice grows: either the value of P 0 decreases and ultimately vanishes in the infinite volume limit, or it becomes approximately constant for large volumes. Both scenarios are shown in Figure 4: for J 1 < J crit 1 P 0 decreases as L with b,c being fit parameters to the data with J 1 < J crit 1 , describe the data very well within their statistical uncertainties. This allows us to conclude that indeed the Polyakov line is a good order parameters for the phase transition in J 1 as it behaves differently on the different sides of J crit 1 , The situation with the Polyakov plane P 01 is more complicated, as it depends non-trivially on both J 1 and J 2 . Moreover this dependence is different in different topological sectors. We show the data in Figure 5. Again, the left panel contains results for J 2 < J crit 2 while the right panel for J 2 > J crit 2 . As opposed to the situation with Polyakov lines, now the plots are no longer similar and there is a nontrivial dependence on J 2 . On the left panel, i.e. for small J 2 , all data sets show a L − 1 2 3 dependence signaling that P 01 vanishes in this region of phase space in the infinite volume limit. This happens no matter what value of J 1 we chose and in both, trivial and non-trivial topological sectors. The right panel contains data for J 2 > J crit 2 . Only a single data set, the blue one corresponding to J 1 < J crit 1 in the topologically charged sector Q 0 = −1, vanishes. In all remaining cases the data show a rather constant value as L 3 is increased, suggesting a non-zero value in the infinite volume limit. Looking from another perspective, in the trivial topological sector Q 0 = 1, P 01 depends only on J 2 , it vanishes for J 2 < J crit 2 and is nonzero for J 2 > J crit 2 , irrespective of J 1 . In the non-trivial topological sector, P 01 vanishes in three corners of the phase space, except of the region where both J 1 and J 2 are large, i.e. J 1 > J crit 1 and J 2 > J crit 2 . Hence, P 01 at Q 0 = −1 is sensitive to both J 1 and J 2 and provides an order parameter for both phase transitions.
Summarizing, for P 01 we have in the limit L 3 → ∞: P 01 > 0 for J 1 < J crit 1 , J 2 > J crit 2 , Q 0 = 1, (3.11) Distinction between phases is seen also by comparing values of P 01 for different coupling constants at one finite value of L 3 , see Table 1.

Summary and conclusions
We have presented an explicit construction of a dynamical lattice model with a local symmetry based on a 2-group. It depends on two coupling constants J 1 , J 2 . We have analyzed the parameter space, first by using dualities to known simpler models, second by simulating the model numerically through Monte Carlo method. Theoretical discussion allows to designate four possible phases in the four corners of the coupling constant plane. In order to study the phase diagram quantitatively, we proposed several candidates for order parameters. Two proposals based on local observables, the average plaquette F and the average cube G, are sensitive to the phase transition only in one of the coupling constants. It follows from the factorization theorem, which we formulate and prove, that F constructed from link variables shows the phase transition in J 1 , whereas G built out of faces shows the phase transition in J 2 . The other two candidates for order parameters are non-local observables. Polyakov lines, which are products of link variables, again, feel only the phase transition in the J 1 coupling constant. Finally, the Polyakov plane exhibits a non-trivial dependence on both J 1 and J 2 and hence can be used as an order parameter for both phase transitions. Furthermore, its expectation value depends on the topological charge sector.
We would like to close this work by mentioning three problems for future study. Firstly, different techniques are required to perform averaging with respect to topological charge sectors. This is because Monte Carlo simulations performed in a fixed topological charge sector do not provide values of weights (partition functions) of distinct sectors. This difficulty is relevant only for those observables for which the average obtained in different topological charge sectors do not agree. The only observable with this property studied in this work is the Polyakov plane. Secondly, it would be interesting to obtain some results about extended surface observables on lattices of topology different than torus, perhaps also for more general crossed modules. Another intriguing question is whether there exists some natural construction of a dynamical higher gauge theory in which factorization theorem does not hold.

A Non-spherical Wilson surfaces
In this appendix we use terminology and notations from [27]. Thus in contrast to the remainder of the paper, this part is not fully self-contained.
It is interesting to ask whether ϕ σ depends on the choice of σ representing the cycle h 2 (σ). If σ is another representative of the same cycle, then σ = σσ 0 for some σ 0 ∈ ker(h 2 ). Thus ϕ σ = ϕ σ ϕ σ 0 . We have to describe ϕ σ 0 . By the characterization of ker(h 2 ) given earlier we have that σ 0 is the product n i=1 (γ i τ i )τ −1 i for some γ i ∈ π 1 (X 1 , * ) and τ i ∈ π 2 (X, X 1 , * ). Thus This element is trivial if either of the following two conditions is satisfied: • ϕ τ i are in ker(∆), i.e. ∂τ i are trivial, • γ i are elements of E which act trivially on Φ; if im(∆) acts trivially (which is satisfied in the model discussed in this paper), this is automatically satisfied if is trivial.
In the language used in the main text, these two conditions correspond to J 1 = ∞ and trivial topological charge, respectively. Assuming that one of these conditions holds, we find that ϕ σ depends on σ only through the corresponding homology class in H 2 (X) (respectively H 2 (X 2 ) if we do not assume flatness of ϕ).