Coarse graining methods for spin net and spin foam models

We undertake first steps in making a class of discrete models of quantum gravity, spin foams, accessible to a large scale analysis by numerical and computational methods. In particular, we apply Migdal-Kadanoff and Tensor Network Renormalization schemes to spin net and spin foam models based on finite Abelian groups and introduce `cutoff models' to probe the fate of gauge symmetries under various such approximated renormalization group flows. For the Tensor Network Renormalization analysis, a new Gauss constraint preserving algorithm is introduced to improve numerical stability and aid physical interpretation. We also describe the fixed point structure and establish an equivalence of certain models.


Introduction
Spin foam models aim at providing a description of the microscopic structure of spacetime and thus a theory of quantum gravity [1][2][3][4][5][6][7][8][9][10][11]. These models can be understood as a non-perturbative definition of the path integral for quantum gravity. To make these path integrals well defined one has to introduce a regularization based on a choice of discretization, i.e. a lattice or, more generally, a triangulation or two-complex. Indeed, spin foams can be understood as generalized lattice gauge theories.
This discrete (a priory auxiliary) structure should not be confused with another feature of spin foams, which is often termed 'Planck scale discreteness' [12][13][14][15][16], namely that the spectra of (kinematical) geometrical quantum observables, like areas and volumes, are discrete. There are thus two different kinds of UV cutoffs, whose interplay has not been fully understood yet. This has to be kept in mind when discussing a possible breaking of (global) Lorentz or (local) diffeomorphism symmetry. A (naive) lattice regularization will generically break these symmetries, see for instance [17][18][19][20][21][22][23][24] for a discussion of these issues in gravity.
We may, however, consider a continuum limit with respect to this auxiliary discretization scale, for example by a coarse graining or blocking procedure, see [25][26][27] for recent examples involving gravity or related to it. A crucial question then is whether Lorentz or diffeomorphism symmetry will be restored in this limit, despite the possibility of still having the second kind of UV cutoff, provided by the discreteness of the spectra, in the theory. That this cutoff does not necessarily lead to a violation of Lorentz symmetry has, for instance, been argued in [28] on kinematical grounds. A full dynamical scenario for 4D gravity where a restoration has been shown to occur is, however, missing, nonetheless see [29,30] for progress in this direction. 1 These questions motivate us to consider a continuum limit which involves many building blocks or large lattices (with many vertices), as this is the limit where one can hope to obtain a diffeomorphism invariant theory. An alternative is the semi-classical limit [36][37][38], in which rather the Planck constant, leading to 'Planck scale discreteness', is taken to zero. To distinguish these two kind of limits we will sometimes refer to the first one as statistical limit.
Experience with other quantum gravity models, such as (causal) dynamical triangulations, has shown [39][40][41][42][43][44][45] that even before the question of restoration of symmetries can be addressed, it is not at all obvious whether such a statistical limit leads to any viable model of spacetime, i.e. whether such a limit will result in smooth four-dimensional spacetime manifolds. Indeed, in this kind of limit statistical considerations become important and it can easily happen that state sums become entropically dominated by configurations not resembling any four-dimensional manifold at all. As we will see, a related issue arises for spin foam models (or other models based on first order/tetrad formulations) where geometrically degenerate configurations might be dominant. Such configurations also turn up in a semi-classical or classical phase space analysis, even if this involves only a single simplex [46][47][48].
Hence it is crucial to investigate which kind of large scale physics or, in other words, phases, are encoded in the candidate quantum gravity models. Phases are often characterized by symmetries, that is such a study might also answer which kind of symmetries might be restored in a large scale / statistical limit. Indeed, making progress in this direction is one of the most pressing issues for the spin foam approach. However, it is also a long standing open issue [49][50][51] charged with a number of conceptual and technical challenges.
One important challenge is the complex structure of the models which lead to very complicated amplitudes as compared, for instance, to QCD. Here, our strategy [11] is to develop a wide range of simplified models which capture essential features of spin foams while being much easier to handle. These simplifications are obtained on the one hand by replacing Lie groups, on which the gravitational spin foams are based, with finite groups. On the other hand, we can also consider 'dimensionally reduced' models (spin net models). In fact, 2D spin net models of the simplest class share many statistical properties with their counterpart 4D spin foam models.
Similar simplified models have been successfully studied, e.g. [52][53][54], to get insights into the large scale behavior of lattice gauge theories. Here we hope for a similar improved understanding of the possible phases that can occur in quantum gravity models. In particular, we will see in the course of the paper, how conjectures or even conclusions for Lie groups can be made based on findings for finite groups. These simplified models can also be interesting in their own right [55][56][57], in particular if an example is found in which some analogue to diffeomorphism symmetry is restored. Indeed topological phases and string net condensates [58] which are studied in condensed matter, also regarding the question of symmetry restoration, are tightly related to 3D spin foams with finite groups [11]. The development of coarse graining and renormalization techniques seems to be the most promising avenue to study the large scale behavior and simplified models allow us to adapt and further refine methods from lattice gauge theory and condensed matter systems. In this work we will therefore apply the Migdal-Kadanoff scheme [59,60] and the tensor network renormalization (TNR) method [61,62]. These schemes involve a regular lattice and, due to this regular structure, they are amenable to efficient numerical simulations.
With this approach we are able to explicitely answer the question of BF symmetry restoration for a range of models and to gain insights into how these results are related to the case of infinite groups. These results should be understood as a first step towards harnessing the power of numerical methods from statistical physics to deepen the understanding of the large-scale physics of spin foam models.
In this work, we will first introduce spin foam and spin net models and write them in ways suitable for coarse graining (section 2). We will also define a particularly important class of models, termed 'Abelian cutoff models', discuss the role of BF / translation symmetry and detail the relationsship between spin foams and nets. Subsequently, in section 3, we discuss the conceptual challenges of coarse graining in this context and argue for the approach pursued here, in particular for the use of a regular lattice.
We then apply the Migdal-Kadanoff and tensor network approximation schemes to coarse grain our models (sections 4 and 5, resp.). In each case, we first introduce the method and highlight some of its analytical properties before presenting numerical results that focus on the question whether renormalization of Abelian cutoff models will restore BF symmetry. In particular, we feature a Gauß constraint preserving TNR algorithm tailored to the geometric interpretation of our models and establish an equivalence amongst certain cutoff models under the TNR renormalization scheme. We conclude by comparing both approximation schemes and pointing out possible future directions of research.

Spin foams and spin net models
Spin foams are a particular class of lattice gauge models (see e.g. [63] for a recent review and [11] for a review emphasizing the relation to lattice gauge and statistical physics models). Such models are specified by variables, taking values in some group G, associated to the edges of a lattice (or more generally an oriented 2-complex) and weights associated to the plaquettes. They can thus also be termed plaquette models.
A related class of models, which will be introduced below, are so called edge or spin net models [11]. Here group variables are associated to the vertices of a lattice (or more generically an oriented graph or 1-complex) and weights to the edges. This class includes the well-known Ising models, based on the group Z 2 . Indeed it will turn out that the structures involved in a spin net model are very similar to those involved in spin foam models -just that where, for instance, weights are associated to 2D plaquettes for spin foams, weights are associated to 1D edges in spin nets, similarly for the group variables and so on. In this sense spin nets are a simpler or dimensionally reduced form of spin foams. There is however one essential difference, which is that spin foams enjoy a local gauge symmetry whereas spin nets only feature a global symmetry, in both cases given by the group G the models are based on. In section 2.3 we will also comment on another relationship between spin foams and spin nets: spin nets can be seen as measuring non gauge-invariant observables in spin foam models.
Spin foams and spin nets are defined by partition functions, and we will first consider a representation of these partition functions as sums over group variables. Via a group Fourier transform we can rewrite these partition functions as sums over variables labeling the irreducible representations of the group G. This is where the name 'spin foam' stems from, as 'spin' refers to the representation labels for the group SU (2). This alternative representation is well known as a duality transformation for both edge and plaquette models, and is usually employed for the high temperature or strong coupling expansion [64]. The models in this representation are not only specified by the dual weights but also by an intertwining projector acting on a certain representation space. Non-trivial spin foam or spin net models can be constructed by choosing this projector to be different from its standard form (which is the Haar intertwiner introduced below) in plaquette and edge models respectively. In the case of Abelian groups we will explicitly construct non-trivial models, the so-called Abelian cutoff models, which only consider representation labels of the group up to a certain cutoff. For these models, we will see that the choice of a non-trivial projector is equivalent to retaining the Haar intertwiner while restricting the dual weights in a particular way. These are precisely the models that we will numerically analyze in sections 4 and 5.
For the non-trivial models, one can then apply the inverse group Fourier transform and again obtain a partition function in terms of group variables. As will be explained below, this will however require the introduction of several group variables per vertex (for edge models) or edge (for plaquette models) [65][66][67]. This representation is termed holonomy representation, both for spin foams and spin nets.
In the next two subsections we will give a short introduction into the main concepts and different representations of spin foams and spin nets. Furthermore we will detail the different possibilities of rewriting these models into tensor network form, as this will be the basis of one of our coarse graining methods, to be discussed later on.

Spin net models
To construct spin net models we start with state sum models formulated over a finite group G and on a graph (one-dimensional complex) with oriented edges. More precisely we consider partition functions of the type where v is the number of vertices in the graph, s(e) denotes the source vertex (starting point) of the edge e, t(e) denotes its target vertex (final point), and the curly brackets under the sum symbol denote that there is a sum per vertex: {gv} = v gv . Here group elements are associated to vertices and weights, which determine the couplings, to edges. Therefore these models are also known as edge models and include the standard Ising model for which G is equal to Z 2 . The weights w e (g s(e) g −1 t(e) ) can be arbitrary functions 2 over the group G. However, if w e are class functions, i.e. invariant under conjugation (w e (g) = w e (hgh −1 ) ∀g, h ∈ G), the model will feature a global symmetry given by the group G: the partition function remains invariant when applying the same conjugation to group variables at each vertex.
The form (2.1) defines the simplest form of spin net models in the representation based on group variables. This representation will be called holonomy representation in analogy to Figure 1: On the left: the three objects associated to every edge. On the right: their schematic representation. Every straight line joining two objects means a contraction of indices.
spin foam models (where group variables are associated to edges and represented holonomies of a connection). Holonomy representations of more general models will require several group variables associated to each vertex, as we will see later.
Via the group Fourier transform, we can change from the above representation in terms of group elements, to the spin net representation, which is in terms of the irreducible representations of the group 3 . Every function on the group can be decomposed in terms of matrix elements of the irreducible representations ρ, being ρ * the dual of ρ. For class functions this decomposition reduces to the character decomposition where χ ρ (g) = dimρ a=1 ρ(g) aa denotes the character. We note that our convention for the delta function over the group, δ G , is Using the property ρ(g −1 ) ab = ρ * (g) ba . (2.5) where we sum over repeated indices. 4 Note that, associated to every edge, there is a coefficient (w ρe ) aebe , and two group representations, ρ e (g s(e) ) aece living on the source vertex and ρ * e (g t(e) ) bece living in the target vertex. The indices of these three objects are contracted, as described schematically in figure 1 .., ρ e4 ) Figure 2: Four-valent vertex with two outgoing edges, e 1 and e 2 , and with two incoming edges, e 3 and e 4 . On the left: representations meeting in the vertex. On the right: schematic representation of the resulting vertex weight.
The first and third groups of indices involve all the edges for which v is the source vertex. The second and fourth groups of indices involve all the edges for which v is the target vertex. Using the schematic representation employed in figure 1, we can represent the vertex weight as in figure 2. Note thatP v can be seen as an intertwiner map, called the Haar intertwiner, acting on a certain representation space for the group G. This representation space, H v , associated to the vertex v, is given by the tensor product of the representations ρ e associated to the outgoing edges and the representations ρ * e associated to the incoming edges The intertwining property of this map is guaranteed by the averaging over the group in (2.8). Indeed, the Haar intertwiner defines an orthogonal projector onto the subspace H inv v of H v , invariant under the group action defined on this representation space.
The same intertwining map will appear in spin foam models. There, the choice of the Haar intertwiner as a projector and face weights to be trivial defines topological models, known as BF -theories. A gauge symmetry, known as translation symmetry, arises in this case, forcing the model not to have local physical degrees of freedom. The analog situation happens in spin net models. The choice of edge weightsw e ≡ id and of the Haar intertwiner as a projector corresponds to the model at zero temperature, with no local degrees of freedom. For the weights w e this amounts to choosing w e ∼ δ G . In this case the projectorsP v are just contracted along the edges of the graph.
Non-trivial models are more interesting. They are constructed by restricting the projector further, i.e. by selecting a subspace of the invariant subspace of H v and by replacing the Haar intertwiner (2.8) by a projector onto this subspace. We will proceed in that way here, and in general assume thatP v is a projector onto some subspace of H inv v . This allows us to obtain interesting models even if we choose the edge weightsw e to be trivial. Indeed, in these models, the original gauge symmetry is broken and local physical degrees of freedom arise. We will see this behavior for the Abelian cutoff models described below.
In general, in order to re-express the partition function Z of such non-trivial models in the holonomy representation, namely as a sum over group variables, we will need to associate more than one group elements to each of the vertices in the lattice. Let us denote the number of edges attached to the vertex v by n v . Now, we assign one group element g(v, e) to each of the edges attached to v and define (2.10) In terms of this vertex amplitude, the partition function reads In case thatP v is given by (2.8), i.e. by the Haar intertwiner, we obtain the form of the partition function given in (2.1), as P v then enforces equality between the group elements associated to one and the same vertex. A different simplification occurs when w e (h) ∼ δ G (h), a choice that we already mentioned above. Then the two group elements g (v,e) and g (v ,e) associated to any edge e have to be equal. Hence, the partition function reduces to a sum over group elements g e = g (v,e) = g (v ,e) associated to edges. This type of models is known as vertex models -the energy of a configuration is now determined by the vertex weights P v ({g e } e⊃v ).
In this work we will apply coarse graining to models with Abelian groups Z q . In this case, all irreducible unitary representations, which we will label by k ∈ Z q , are one-dimensional and defined by their characters χ k (g) = exp( 2πi q k · g) for g ∈ Z q . The transformation between functions on the group Z q and on the dual, equal to the space of characters, which is also given by Z q , is given by the discrete Fourier transform (2.12) Characters for Abelian groups are multiplicative, i.e. χ k (g 1 · g 2 ) = χ k (g 1 ) · χ k (g 2 ), and also χ k (g −1 ) = χ −1 k (g) = χ k (g). Moreover, the delta over the group is now the q-periodic delta. It verifies q −1 g δ (q) (g)f (g) = f (0). The spin net representation of the partition function simplifies in the case of Abelian groups where ε e v is equal to +1 (−1) if v is the source (target) of e. The projectorP v implements the Gauß constraints at the vertex v. It indeed projects to the irreducible subrepresentation in the tensor product of all representations associated to the outgoing edges and the tensor product of dual representations of incoming edges. The tensor product is one-dimensional and equal to the trivial representation if the oriented sum of the representation labels k e is equal to zero.
Note that the spin net representation is the starting point for the high temperature expansion [69]. The infinite or high temperature fixed point is represented byw e (k) = δ (q) (k). On the other hand the zero or low temperature fixed point is given byw e (k) ≡ 1. This corresponds to weights w e (h) = δ G (h) in the original group representation. The model is 'frozen', i.e. the group elements at the different vertices have to agree (assuming that the graph has only one connected component).
As commented before, for this choice of weights there is a gauge symmetry. It is associated to the faces, i.e. the two-dimensional cells (here we are assuming that the graph is actually given by an orientable 2-complex, i.e. the 2-dimensional cells are well defined). Associating to every face f an element k f ∈ Z q we define a gauge transformation acting as where ε f e is +1 (−1) if the orientations of e and f agree (disagree). Under such a transformation the contribution of a configuration {k e } e to the partition function Z does not change. Choosing either the edge weightsw e or the vertex projectorP v to be non-trivial will in general break this translation symmetry either completely or down to a smaller symmetry. Choosing a non-trivial projectorP v will in general result in vertex models, as one can basically reduce the set of vertices allowed by the Gauß constraints even further.
We are now in position to introduce the Abelian cutoff models. As said before, for the low temperature fixed point (analog to BF in spin foam models) the weights are The Abelian cutoff models are derived from this by 'cutting off' the sum at some value K such that some of the dual weightsw e vanish. Explicitly, where we will consider only even q, hence K ≤ q 2 . Here the range for k is given by −q 2 < k ≤ q 2 . Also, note that the symmetry conditionw(k) =w(−k) ∀k is fulfilled. This requirement is desired in the quantum gravity setting because it ensures that the model does not depend on edge orientations and certain types of face and edge subdivisions [70,71].
Abelian cutoff models could be equivalently described by a restriction of the projectorP v (and keeping the weightsw e (k) ≡ 1). They provide a simple example of breaking the translation gauge symmetry from the frozen model in order to introduce physical degrees of freedom. Also this choice of model corresponds to a regularization one often choses for BF lattice theories with Lie groups [72]. In this case the orbits of the translation gauge symmetry are non-compact and the evaluation of the partition function gives generically infinity. Different methods of regularization have been developed, one would be equivalent to introducing a cutoff K (for a theory with group G = U (1)). One can then ask whether these regularized models would flow back to the full BF model, or more generally the low temperature fixed point, under coarse graining. We will consider this question in sections 4 and 5.
So far we have presented the partition functions for spin net models as sums over group variables (the holonomy representation) or representation labels (the spin net representation). An alternative is to write the partition function as a contraction over tensors attached to the vertices of the underlying lattice (or some graph associated to the lattice), that is, in the tensor network representation. The tensor network representation is commonly employed in statistical and quantum systems [61,[73][74][75][76], since it is especially suitable for developing techniques of renormalization. We will make use of this representation of the partition function in section 5, where we will apply the renormalization approach to spin net models with Abelian groups.
In this representation, the contraction of the indices is prescribed by the edges. The partition function can hence be expressed as a tensor trace   This way of representing the partition function is related to the so-called graphical calculus [77][78][79][80], which is often employed in the spin foam literature.
In the particular example of Abelian models in the spin net representation, we first absorb the edge weightsw ke into the vertex weightsP v ({k e } e⊃v ) by distributingw 1/2 ke factors to each of the adjacent vertices, namely in every vertex we define the tensor The partition function is then given by the tensor trace (2.17) (see figure 3 and figure 4). More generally, we can find tensor network representations also for the non-Abelian models in the different representations. The holonomy representations is a convenient starting point for a low temperature expansion, the spin net representation to the high temperature region. For the spin net representation we can proceed as for the Abelian models. For the holonomy representation assume that we have edge weights such that w e (h) = δ G (h). We then just need to understand the group elements g e as indices attached to the tensors T s(e) = P s(e) and T t(e) = P t(e) to see that (2.11) can be rewritten as a tensor trace. For the more general case of non-trivial edge weights we can introduce another set of rank two tensors w to the midpoints of the edges.

Spin foam models
Spin foams are in many respects similar to spin nets. The main difference is that spin foams are gauge theories, formulated with a gauge group G, which here will be a finite group. Furthermore spin foams require an oriented two-complex for their definition. This implies, that we have a well defined notion of oriented edges as well as oriented faces (which are the 2D cells of the complex, or the plaquettes for a regular lattice).
In order to define a spin foam model, we assign a group element g e to every edge e of the two-complex and a weight w f : G → C to every face f . The state sum model is then defined by the partition function (2.19) where e denotes the total number of edges in the two-complex.
The function w f is a class function. Furthermore, in (2.19) w f depends on the group elements only through the holonomy h f around the closed loop of edges forming the face f , that we will call curvature. Let us make its definition explicit. For that, we recall that the relative orientation between a face f and any edge e in the boundary of f is denoted by ε f e , and it is equal to +1 (−1) when e and f have the same (opposite) orientation. Given a face f bounded by the ordered sequence of edges e 1 , e 2 , · · · , e n , the associated curvature is given by as depicted in figure 5a. These properties of the weight w f guarantee that the partition function is invariant under gauge transformations g e → g s(e) g e g −1 t(e) . As before, s(e) denotes the source vertex of the edge e while t(e) denotes its target vertex.
The partition function (2.19) describes standard lattice gauge theories. The weights can for instance be chosen to emulate the Wilson action where α is a coupling constant. For the choice w f (h f ) = δ G (h f ), with δ G the delta function over the group defined as before, the partition function only sums over (locally) flat holonomies.
This is a discretization of BF theory, which is a topological field theory (without propagating degrees of freedom). It coincides with the zero temperature fixed point or zero coupling fixed point of lattice gauge theories of Yang Mills type. In order to obtain a representation of the partition function (2.19) as a sum over representation labels (spin foam representation), we again apply the group Fourier transform. Here we only need to decompose class functions into characters (2.23) We introduce this decomposition in equation (2.19) and individually carry out the sums over the group elements (note that given an edge e, the groups elements g e and g −1 e appear in f w f as many times as number of faces share the edge e). We obtain the following expression for the partition function where, associated to every edge, we have defined the projector In the above tensor, the first and third groups of indices (distinguished by the superindex f + ) involve all the faces that have e in their boundary with the same orientation as the face. These groups of indices have therefore as many indices as number of faces with ε f e = 1. In turn, the second and fourth groups of indices (distinguished by the superindex f − ) have as many indices as number of faces f with ε f e = −1. We show an example in figure 5b. To every pair face-edge two indices are associated, a f s and a f t , that we attach to the vertices of f that bound the edge e, namely s(e) and t(e). In equation (2.24), for every face f of the two-complex, there is a sum for every vertex v belonging to that face. Note that every index a f v appears twice, since there are two edges meeting at the vertex v and bounding the face f . Then, the product over edge projectors contracts all the indices.
Associated to every edge there is a Hilbert space H e . Let us denote by f 1 , · · · f m the faces attached to e. Then, Here, to keep the notation simple, we are assuming that ε f e = 1 for all the faces 5 . The corresponding tensorP e defines an orthogonal projector onto the gauge-invariant subspace of the edge-Hilbert space H e (P e ψ) a 1 a 2 ···am := (P e ) a 1 a 2 ···am; for ψ ∈ H e . Here, as for the spin net models, the projector (2.25) is the Haar-intertwiner. As before, more general spin foam models are constructed by restricting the Haar intertwiner to proper subspaces of the gauge invariant subspace of H e . For spin foam models one usually reorganizes the partition function (2.24) such that amplitudes can be associated to vertices. To this end, one decomposed the projectorsP e by introducing an orthonormal basis ι e k , k = 1, . . . , m, for the invariant subspace of H e . By adjusting the basis we can decompose any gauge invariant projector as (here m ≤ m and m = m for the Haar intertwiner) so that the indices attached to t(e) are now carried by the intertwiners |ι e and the indices attached to s(e) are carried by the intertwiners ι e |. According to the description above, we can contract within every vertex the corresponding ι e obtaining vertex amlitudes A v (ρ f , ι e ), that depend on the representations ρ f and intertwiners ι e associated to the faces and edges that meet at v. With this, the partition function can be written in terms of vertex amplitudes via This is the usual description employed in spin foam models for quantum gravity.
As for the spin net models, we can also define a holonomy representation [65][66][67] for spin foams using the inverse group Fourier transform. In case we are dealing with a non-trivial edge projector (2.28), i.e. not the Haar intertwiner, the resulting holonomy representation will be of the form of a generalized lattice gauge theory. That is, instead of only one group variable associated to every edge we will have as many group variables attached to a given edge as there are faces attached to this edge. Furthermore, there will not only be face weights w f but also edge weights P e which result from the transformed edge projectorsP e .
Let us now consider the case of Abelian groups Z q . As for spin net models, the spin foam representation simplifies since the irreducible representations are just one-dimensional and the group Fourier transform is just given by the usual discrete Fourier transform. Namely, we obtain where the q−periodic delta function δ (q) (·) enforces the Gauß constraints, now based on the edges (instead on the vertices as for spin nets). In the last step we just splitted the delta functions over every edge into two delta functions over every vertex, in order to define the vertex amplitudes (which here are A v = e⊃v δ (q) f ⊃e ε f e k f ). In case the face weights are given byw k f = 1, which is the BF theory case, we obtain an additional symmetry for the partition function (2.30), which as commented before is known as translation symmetry. For spin foams the gauge parameters k c ∈ Z q are based on the 3-cells of the lattice (here we assume that the 3-cells are well defined, for a regular hypercubic lattice these would be the 3D cubes): where ε c f = +1 (= −1) if the orientations of the 3-cell c agrees (disagrees) with the one of the 2-cell f . As for the spin net models this translation symmetry does not change the validity of the Gauß constraints appearing in the partition function (2.30). For the gravitational spin foam models (in 3D) this type of gauge symmetry gives the diffeomorphism symmetry underlying general relativiy [20,31]. It hence has a special status and is indeed deeply intertwined with triangulation or more generally discretization independence of the models [20-22, 25, 27].
For coarse graining we will consider the Abelian cutoff models with face weights coinciding with the edge weights in spin net models (2.16). For the cutoff models the translation symmetry will be broken as the weights are now no longer constant in k f . The motivating question will be to see whether these models flow back to the BF phase, in which the translation symmetry is restored.
Finally let us note that also spin foams can be represented in various ways as tensor networks. One possibility would be to start with the representation (2.29) involving vertex amplitudes. To obtain an algebraically similar form to the tensor network representation of spin net models one would however keep the edge projectorsP e intact. To this end we associate the tensorT e =P e to the midpoints of the edges of the lattice. The edges of the lattice carry a number of indices which are all contracted with each other in the lattice vertices according to the description below (2.25). Furthermore we have to take care of the face weightsw ρ f and the sum over representation labels ρ. This can be achieved by introducing another type of tensorT f =w, see figure 6. If we work with a cubic lattice these tensors are four-valent carrying as indices representation labels: The second factor is equal to one if all representation labels in the argument coincide and zero otherwise. This tensor is connected by auxiliary edges to the four adjacent edge tensors T e ensuring that the sum over the representation labels involves the same representation label for every face.

Relation between spin foam and spin net models
Here we want to comment on an interesting relationship between spin net and spin foam models. Namely, one can understand spin nets as expectation values of (gauge symmetry breaking) observables inserted into the spin foam partition functions. We will illustrate this only for the simplest case: the spin foam represents BF theory and the spin net is of the form (2.1), i.e. the vertex projector is given by the Haar intertwiner.
Indeed, the partition function (2.1) can be rewritten in the form where f runs over a set of faces whose boundary edges and vertices generate the initial graph over which the spin net model is defined. Moreover, the two-complex made up by this set of faces must be simply connected for equation (2.32) to be valid (otherwise we have to amend the condition that holonomies along non-contractible loops should be trivial). The constant of propotionality comes from the normalization of the delta functions and can be absorbed by a redefinition of the weights w e (g e ) in the second line. The second line of equation (2.32) is the result of introducing the product of edge weights w e (g e ) in the partition function for BF theory (w f = δ G ) on a simply connected two-complex. Similar observables have been considered in the context of 3D quantum gravity, namely for the Ponzano-Regge model [81,82] and have been interpreted as Feynman diagram evaluations. That is, the edge-weights w e (g e ) can be understood as propagators and the spin net model as a Feynman diagram evaluation.

Coarse graining methods
Having introduced the models of interest in the previous section, we now turn to the challenge of implementing a coarse graining procedure. In this section, we will outline our general approach and discuss the conceptual issues that arise. The application of two approximation schemes in our context will then be discussed in the following sections.
Although gravitational spin foams are usually defined on a general triangulation or twocomplex we will here consider coarse graining on a regular lattice. This allows us to actually make explicit computations and to use methods from lattice gauge theory and condensed matter systems. One might object that using a regular lattice will introduce a background structure, spoiling background independence. There are, however, indications [25,27,83] that a restoration of diffeomorphism symmetry will be connected with a notion of triangulation or discretization independence, and hence the choice of a particular underlying lattice may not matter. Indeed, the universality phenomenon of statistical systems also suggests that the details of the chosen lattice might not matter for the questions we are interested in, i.e. a characterization of the possible phases of the models, or whether spin foams can avoid degenerate phases. Nevertheless one should study whether the results depend on the choice of lattice.
The development of a scheme where order parameters or coupling strengths might be locally varying is another conceptual challenge [49,50], which we will not address here. This would be appropriate for a random lattice or for situations with a very inhomogeneous dynamics. Again, we think that developing feasible coarse graining methods for spin foams and nets on a regular lattice is an indispensable first step.
Furthermore, for gravitational systems, a regular underlying lattice can nevertheless represent very inhomogeneous or irregular geometries. This is due to the fact that the dynamical variables are the geometrical data which also determine the lattice geometry and moreover the physical scale. In this sense a very fine lattice can carry very different (coarser) geometries [20,84]. 6 This opens up the possibility that refining in lattice quantum gravity is even equivalent to summing over (a class of coarser) lattices [84], thus eventually leading to a discretization independent theory at the fixed point [25,27,83]. Indeed, systems with diffeomorphism (like) symmetry might add interesting new insights to the theory of phase transitions.
Another issue, which has to be addressed for any coarse graining or renormalization scheme, is which space of models one considers the renormalization flow to take place in [51]. Indeed, this question is not quite obvious for spin foam or spin net models as, for instance, spin net models mix aspects of edge models, where couplings are along edges, with vertex models, where couplings are based at the vertices. type of degrees of freedoms couple in a certain way might not be available for the non-trivial models.
Furthermore, spin foam models are constructed to be as independent as possible from the underlying discretization. This translates into certain invariance properties of the amplitudes under a certain class of subdivisions, those which effect edges or plaquettes, but do not lead to an increase in the number of (non-trivial) vertices [70,71,86,87]. Spin foam models, which are constructed using trivial face weights but non-trivial projectors, will satisfy this invariance. Invariance under this subclass of subdivisions can be seen as a first step towards a discretization independent model. One might ask whether it is possible to come up with a renormalization scheme in which this form of the spin foam amplitudes and the invariance properties under face and edge subdivisions is preserved, see also [67]. However, this seems not to be very likely, as long as the amplitudes are not invariant under all kinds of subdivisions and hence renormalization flow is trivial or at a fixed point. An intuitive reason is that the faces and edges of a coarse grained spin foam are 'effective' building blocks, containing a huge number of bare vertices, edges or plaquettes. That is, a change of the effective triangulation, even if it only involves a subdivision of edges and faces, will correspond to a more complicated change of the underlying 'bare' triangulation. Indeed, we will find that under the renormalization schemes presented here, the invariance property under edge and face subdivisions is not preserved.
A general problem with real space renormalization schemes of (higher than one-dimensional) statistical models is that an increasing number of non-local couplings appear with each blocking step. 7 To keep the renormalization flow in a space with a finite number of parameters, some sort of truncation or approximations scheme has to be used. We will consider two such schemes. The first one, the Migdal-Kadanoff [59,60] scheme, is based on a truncation to local couplings. The derivation of this scheme is based on arguments that rely on the standard form of edge and plaquette models. Here an important question for future work would be to generalize this scheme to proper spin foam or spin net models, which constitute a rather generalized form of these models.
Whereas the Migdal-Kadanoff scheme is based on a blocking of the group variables, the second scheme, based on tensor network renormalization, is blocking the degrees of freedom encoded in the representation labels. This might be an advantage as the representation labels are related to metric degrees of freedom in the geometric spin foam models. Hence, this kind of blocking is much closer to the blocking used in the gravitation models in [25,26], which is derived by geometrical arguments. Moreover we will see that this scheme allows direct access to the behavior of the vertex and edge projectors and makes use of the properties of representation theory encoded in the models.
The tensor network renormalization can, however, also be applied in the holonomy representation of the spin foam and spin net models, see also [90]. Which scheme to use might depend on which region one aims to explore, the high temperature (where the spin net/ spin foam representation is more appropriate) or the low temperature region (where the holonomy 7 Indeed, such non-local couplings are essential to regain diffeomorphism [26] or Lorentz symmetry [88,89]. This at least applies to models with genuine field degrees of freedom, i.e. 4D gravity, but not for 3D gravity, which is a so-called topological field theory. Nevertheless, 3D gravity (which can be described by BF theory, that is the low temperature fixed point of lattice gauge theories) is an important test case for which the restoration of symmetries can be studied in a simpler setting.
representation might be more useful). Indeed, tensor networks are an extremely general tool, in which it is also possible to consider the emergence of different kinds of effective degrees of freedom, long range order as well as topological phases [91]. Another advantage is that tensor networks can handle non-positive weights [90]. This is an important point for gravitational spin foams where oscillating amplitudes appear, preventing the use of Monte Carlo simulations.
The issue of non-local couplings in this renormalization scheme appears in the form of having more and more degrees of freedom. One has to choose a cutoff for this number. The advantage as compared to the Midgdal Kadanoff scheme is that the accuracy of this scheme can be systematically improved by increasing this cutoff. Furthermore, it is possible to study the dependence of the results on the choice of cutoff.
On the other hand the tensor network renormalization scheme requires much more effort than the Migdal-Kadanoff scheme and offers less analytical control. This is of course related to the much bigger parameter space one is considering in the case of tensor networks. A crucial question for future research is how feasible the tensor network scheme will be for higher dimensional systems, as most work performed so far is for one and two-dimensional systems. Our application of the tensor network renormalization method is also restricted to two-dimensional spin net models.

The Migdal-Kadanoff approximation
Migdal-Kadanoff (MK) approximations are a simple tool to overcome a key difficulty in realspace renormalization: the introduction of non-local couplings in the renormalized action (which in the statistical physics language is the Hamiltonian) with each coarse graining step.
When coarse graining the one-dimensional Ising model with nearest-neighbour interactions by decimation, the state sum factorizes. Integrating out the chosen spins gives a renormalized Hamiltonian of the same form as the one we had started with [60]. However, this no longer happens in higher dimensions: generically, the renormalized Hamiltonians feature new, longerranged interaction terms that render it impossible to continue with the procedure [59]. The central idea behind the MK approximation is to substitute this renormalized Hamiltonian by an approximate Hamiltonian which features the same interaction terms as the original Hamiltonian, thus making the renormalization transformation form-invariant.
The original proposal [59] of Migdal amounts to neglecting certain terms in the state sum, while strengthening others. Applied to the two-dimensional Ising model, this translates to moving bonds (couplings along the edges) away from those nodes which are to be eliminated by decimation. In this way the valency of those particular nodes is reduced to two and the situation is thereby effectively reduced to the one-dimensional case (see figure 7).
Kadanoff subsequently identified bond-moving as a special case of a general approximation scheme for which he derived a bound on the free energy [60]. Further justification for this approximation is derived from Monte Carlo simulations [52,53]. In practice, MK approximations do fairly well at finding fixed points and phase transitions and not so well at determining the order of these transitions [54]. However, specific fixed points of Kosterlitz-Thouless type might not be reproduced by this approximation [92]. MK approximations are computationally efficient and comparably easy to implement but not refineable in a systematic way. based on the Abelian group Z q . Bond-moving then involves dropping some of the terms in (4.1) and replacing some of the others with

Migdal-Kadanoff relations explored
leading to convolution of Fourier coefficients, the main feature of MK approximations. The subset of group variables corresponding to dropped terms can then be integrated out without generating non-local couplings. One can proceed similarly for lattice gauge theories or plaquette models of the form (2.19). For both kind of models different versions of bond-moving and different decimation schemes can be defined. In particular, in anisotropic methods the renormalized couplings will be different in different lattice directions, while in simpler isotropic methods the couplings will remain homogeneous (see figure 7). We will only consider the isotropic methods here, which can also be made exact on so-called hierarchical lattices [93].
For the normalized weights Q(k) :=w(k)/w(0) this results in general recursion relations of the form [92] Here, γ and λ are constants specific to the model and derivation in question. However, models with fixed γ · λ are qualitatively equivalent and share the same fixed point structure in the following sense: Two instances of (4.3), with exponents γ, λ;γ,λ resp. and γ · λ =γ ·λ are related by given initial configurations that satisfỹ Typical values for the exponents are γ = 1 and λ = 2 for 'isotropic' 2D spin net and 4D lattice gauge models and γ = 1 and λ = 4 'isotropic' 3d lattice gauge theories. The corresponding fixed point equations are given by The high and low temperature fixed points ( Q(0) = 1, Q(k > 0) = 0 and Q(k) = 1 ∀k, resp.) are common solutions of eq. (4.6), independently of the exponents. Furthermore, for some fixed factor d of q, there is a class of fixed points given by which derive this property from being invariant under d-fold cyclic permutations. More generally, for a fixed factor d of q, (symmetric and normalized) configurations with Q(k) = 0 for k mod d = 0 and those given by form invariant submanifolds of the parameter space. There are also more nontrivial fixed points, including the higher dimensional analogue of the Ising fixed point on the one-dimensional invariant line connecting the low and high temperature fixed points. For Z 2 , this point is predicted to be a nontrivial solution of In the 2D standard Ising model case (with exponents γ = 1, λ = 2) the solution is given by Q = w 1 /w 0 = 0.296 which corresponds via kT = Artanh(Q) −1 to a temperature of kT = 3.282 (the exact solution is given by kT c = 2/ log(1+ √ 2) ≈ 2.269 [94]). Also note that MK approximations only predict this Ising-type fixed point for exponents with γ · λ > 1. Figure 8 illustrates the afore-mentioned features for 2D spin net and 3D spin foam models with group Z 4 .

Restoration of BF symmetry in Abelian cutoff models
Cutoff models are certain initial configurations for the MK renormalization group flow which are of particular interest to quantum gravity because they are derived from BF theory. The question of interest here is whether the topological nature of BF theory, which is destroyed in the construction of the cutoff models, will be restored under the renormalization group flow. With a possible restoration of the BF phase, also the translation symmetries (2.14)-(2.31) would be restored, which in the 3D gravity models correspond to diffeomorphism symmetry. Whereas the BF phase -or low temperature fixed point (LTF)-in a 3D gravity setting represents flat space time, the high temperature fixed point (HTF), at which Q(0) is equal to unity and vanishes for all other labels k, corresponds to a geometrically degenerate phase, in which all geometric (length) observables have vanishing expectation value. However, whether the models flow back to BF or not depends very much on the initial configuration in question. Instead of focussing on one model with one specific group, we here address the differences between different models for finite Abelian groups Z q with varying q, leaving other (non-Abelian) groups for further research.
Let us first consider the case of 3D lattice gauge theory / spin foam models over Z q for varying q where the exponents in (4.3) are given by γ = 1, λ = 4. Here, different initial configurations parametrized by q and K converge quite fast (after 5 to 10 iterations) either to the HTF or to LTF, see figure 9b. In general, there are two competing effects encoded in the recursion relations (4.3): the convolution leads to a broadening of the function Q(k), whereas the exponents γ and λ lead to a dampening effect (as the Q(k) ≤ 1).
For the 3D gauge models most configurations flow to the HTF, that is the dampening factor is quite strong. Regarding the question of restoration of the translation symmetries we see that this happens only in cases where the initial configuration is already quite close to the BF configuration, which coincides with the LTF.
These observations are in accordance with similar, but analytical work done in the generalized case of U (1) [92]. There it is found that for 3D U (1) gauge models all configurations satisfying certain conditions 8 flow to the HTF.
These results have been extended to 3D lattice gauge models with non-Abelian compact groups U (N ) and SU (N ) by [95]. Hence this applies also to the 3D gravitational (Ponzano-Regge) model, which is based on SU (2) (assuming a change from a lattice based on tetrahedra to a cubic lattice does not matter) and we have to conclude that within the Migdal-Kadanoff approximation translation symmetry / diffeomorphism symmetry is not restored. Note, however, that for finite groups (i.e. finite q), there are some (even) cutoff configurations that do flow back to the LTF or BF phase. Here we might draw the conclusion that it is easier to restore a compact symmetry as compared to a non-compact one, as the translation symmetry is based on a compact parameter space for Z q and on a non-compact one for proper Lie groups. Here it would be interesting to see whether for some analogous modifications of the Tuarev-Viro models [33], which describe 3D quantum gravity with a cosmological constant and are based on quantum groups, such a restoration of the (here compact) translation symmetries occurs. See also a discussion of related issues in loop quantum gravity quantization of 3D gravity with a cosmological constant [34,35] and a classical coarse graining treatment of the same system, where diffeomorphsim symmetry can be restored [25]. To study this question for the Tuarev-Viro models one would have to adjust the Migdal-Kadanoff method to quantum groups or to apply alternative renormalization schemes, such as the tensor network scheme described in the next section.
For 4D 'isotropic' lattice gauge models, or equivalently for 2D spin net or edge models, the exponents in (4.3) are given by γ = 1, λ = 2. Hence the dampening effect is much less pronounced than for the 3D gauge models. Indeed, now most configurations flow back to the LTF or BF phase, see figure 9a. An exception are the K = 1, q ≥ 8 configurations which flow and would be stable at least up to the number of digits displayed in (4.10). This configuration converges after 970 iterations to the low temperature fixed point. 9 Note however that this number of iterations corresponds to an extremely large lattice (in lattice units). This type of behavior is typical for fixed points with unstable directions but could also occur for quasi fixed points. To differentiate between these two cases, we considered a one-parameter deformation of the q = 14, K = 2 model, for which we changed the Q(1) = Q(2) = Q(12) = Q(13) values from unity to an arbitrary parameter 0 < x < 1. Indeed, there is a phase boundary for x ∼ 0.365 which leads to a non-trivial fixed point (we give only the first two non-vanishing digits) This fixed point is considerably different from the configuration (4.10) and also the initial configurations parameterized through x = 1 and x = 0.365 are quite different. Hence we see that a rather large portion of parameter space is dominated by the unstable fixed point. To obtain convergence to either LTF and HTF extremely large iteration numbers are necessary. In other words, the phase transition between LTF and HTF is very weak and the phase boundary not very pronounced.
Indeed, in two-dimensional systems and in the limit of a continuous symmetry group, such as U (1), one cannot expect the usual type of second order phase transition between the symmetry breaking phase (which here would be LTF) and the disordered phase (HTF). This is explained by the Mermin-Wagner (Coleman) theorem stating that continous symmetries cannot be spontaneously broken at finite temperature [96][97][98]. Nevertheless, there are two phases, both in the two-dimensional system with global U (1) symmetry [99,100] connected by a Kosterlitz-Thouless transition [101] and in the four-dimensional gauge system [102,103]. This is a phase transition of infinite order. However, for the (isotropic) MK relations it was proven by Ito (again under certain assumptions on the initial configurations, including positivity of weights both in group and Fourier space), that this phase transition is not detected [92]. But also here extremely slow convergence (of all configurations towards HTF) has to be expected due to the existence of quasi fixed points. Our findings are explained by these considerations (although we use initial configurations not satisfying Ito's assumptions), that is by going to larger q we have to expect weaker and weaker phase transitions. Furthermore, the Ising type fixed point on the invariant line between LTF and HTF is drawn continuously towards HTF, indicating the absence of the transition in the limit of large q. For 4D lattice gauge theories with non-Abelian Lie groups it is speculated [92] that the MK relations provide a better approximations than for Abelian ones. The confinement conjecture states that in contrast to Abelian groups there is no phase transition for systems with non-Abelian Lie groups [103].

The tensor network renormalization scheme
Here we will discuss a second coarse graining method and apply it to spin net models in 2D in the spin net representation.
We have shown in section 2.1 that the spin net models can be easily brought into tensor network form. For these a number of real space renormalization techniques have been developed in recent years [61,62,104]. Moreover, tensor networks have became popular not only as a tool to formulate partition functions for 'classical statitistical models' but also to provide a variational ansatz for trial wave functions [73,76,105] for quantum statistical models. For a variational ansatz one has to find the expectation value of the Hamitonian with respect to the trial wave functions. The computational techniques [105] are similar to the tensor network renormalization group techniques. This is another point that motivates us to consider renormalization of spin net models, as structures very similar to spin nets, the so-called spin networks, appear in the canonical or Hilbert space formulation of spin foam models. Hence the renormalization techniques considered here could also be useful to find e.g. the physical wave function (ground state) via a variational ansatz. Indeed, the tensor network ansatz has also been developed to describe topological phases [62,105], which often are represented by the so-called physical wave functions of BF -theories, which are the starting point of spin foam quantization. Moreover, as we have also seen in section 2.1, tensor networks and spin networks are naturally related [74,78]. More generally, tensor networks are a general tool for graphical calculus [78,79], which becomes especially powerful for representation theoretic models such as spin foams and spin nets [77,80].
The starting point of the tensor network renormalization (TNR) method is to write the partition function of a given model as a contraction of tensors associated to vertices of a graph (or lattice) T abcd T aef g T bhij · · · . (5.1) The contraction of indices is along the edges of the graph, that is every edge of the graph carries one index. As these are summed over we can interpret the indices as the variables or degrees of freedom of the model. The dynamics is encoded in the choice of tensor. The tensor network itself, that is its underlying graph, can also be interpreted as a lattice in space or space-time. Choosing appropriate subsets of tensors and contracting all tensors inside each subset according to the connectivity given by the network will result in a network with fewer vertices and edges. Hence, this procedure corresponds to the blocking procedure in real space renormalization. Each subset results in a new 'effective' tensor T describing an effective model. Notice, however, that in the two and higher dimensional case the effective tensors T generically 10 carry more indices than the original tensors T . For a regular lattice the number of indices associated to the effective tensors T grows with the number of iteration steps. For the underlying graph it means that the valency of the effective vertices will grow -mostly in the form of having multiple edges between pairs of vertices. These multiple edges can be summarized into effective edges, which then carry indices with an exponentially growing range.
Here is where one has to choose an approximation such that the indices run over a prechosen maximum number of values D c . By increasing D c the approximation can be improved systematically. Ideally, this approximation should pick out only the relevant physics and neglect the irrelevant short distance fluctuations. The details of how this selection is implemented depend on the scheme, here we will follow a refined version of [61,62].
The TNR method is very general as many different models can be written in tensor network form. That is in principle one can also flow between different models based on different kind of variables. models we obtain the same initial tensor network models for which also follow the same renormalization flow. On the other hand one loses a direct physical interpretation of the 'blocking procedure', i.e. how the effective/ blocked degrees of freedom are built up from the microscopic ones. This information can be supplemented by studying expectation values of (coarse grained) observables. Below we will introduce a method which allows to keep some physical interpretation of the indices associated to the effective tensors. This is however specific to the spin net/ spin foam representation, where the indices are group representation labels and for the gravitational spin foam models carry geometric information, such as lengths and area values. Hence we can argue that this method would correspond to blocking over these area or lengths variables. Such a blocking procedure was also employed in [25,26] for the classical Regge model and has the advantage of a direct geometrical interpretation of the coarse graining procedure.
Here we will consider the TNR method for a regular lattice but the principle is also applicable to generic graphs arising for instance from random triangulations or Feynamn graphs. In this case one needs, however, to find some suitable approximation scheme to prevent an exponentially growing index range of the effective tensors.

Gauß constraint preserving TNR method
In this section we will shortly describe the TNR method following [61,62] applied to 2D Abelian spin net models. We will, however, introduce a technique to keep the Gauß constraints explicitly valid throughout the renormalization process, see also [74,75]. The reason for doing this is that the Gauß constraints have an immediate geometrical information: In a given spin net (with oriented edges) consider any region such that its boundary cuts only through edges. Then only those configurations will contribute to the partition sum for which the sum of all ingoing indices is equal (modulo q) to the sum of all outgoing indices. This means that the Gauß constraints should also hold at the effective vertices, which arise from blocking all the vertices in certain  regions. We will first review the method for a general 2D tensor network model based on a square lattice and afterwards specify to the case of spin net models and deal with the Gauß constraints.
Consider a 2D tensor network based on a square lattice, so that the tensors T abcd are of rank four, see figure 10a. An obvious way to proceed would be to contract always four tensors along a square and to define in this way a new effective tensor which would now carry four double indices.
However, to find a suitable approximation, that is a method to keep the index range constant, one proceeds differently. The first step is to decompose the tensors T into a product of two other tensors S. This is performed in two different ways according to the partition of vertices into odd and even ones. A vertex is even, respectively odd, if the sum of its lattice coordinates is even, respectively odd.
For even vertices we decompose (see figure 10b) with positive singular values λ i and unitary matrices U and V . We can then define S ab,i If we keep the range of i as in equation (5.3) the index range of the tensors T would grow exponentially with the number of iterations. This is where the key approximation step comes in, namely to consider only the D c largest singular values in the decomposition (5.3). This approximation is justified as the partition function is a trace over the tensors, thus involving the sum over the singular values. The validity of the approximation can be checked by comparing the values of the neglected singular values against the largest singular values in the SVD [61]. One can choose a rescaling after each iteration step such that this largest singular value is equal to one. Implementing the cutoff D c in the number of singular values in the decomposition (5.3) we will obtain a flow in the space of tensors of rank four with a constant index range given by D c .
The SVD does not only serve as an approximation method but leads also to a field redefinition. Here the field variables are given by the indices over which the tensors are contracted. In the SVD these tensors are linearly transformed, which also induces a transformation on the fields. The transformations aim at an efficient representation of the partition sum, i.e. involving a minimal range of indices or equivalently minimal number or range of variables. The SVD is not unique, in particular for degenerate singular values one can add rotation matrices acting on the eigenspaces associated to the degenerate singular values. For the spin net models we will, however, modify the method so that the kind of field redefinitions that can occur are restricted. This is related to preserving the Gauß constraints throughout the renormalization process. It also has the advantage that the indices keep their original physical interpretation, which for the gravitational models are related to geometrical quantities.
Let us now specify to the Abelian spin net models in the spin net representation. Here it is convenient to introduce an orientation for the edges: for the square lattice we will choose all horizontal edges to point to the right and all vertical edges to point upwards. In this case the initial tensor T is of the form (the indices are anti-clockwise cyclically ordered starting from the leg pointing to the right, as in figure 10b) where u(·) = w(·) in the notation of (2.16) and a = 0, . . . , q − 1 for a model based on Z q . The delta function factor signifies the Gauß constraints. Because of this factor the matrices M 1 and M 2 can be brought into block diagonal form, namely for M ab,cd 1 we have the condition that for non-vanishing entries, whereas for M cb,ad Here the indices i, j label the non-vanishing blocks.
In the first iteration step the decomposition into the tensors S can be obtained exactly and involves at most q non-vanishing singular values We assume that D c > q (or alternatively that D c is bigger than the number of non-vanishing u(a)), so that no approximation is necessary at the first iteration step. Note that at least D c = q is necessary to flow to the low temperature fixed point, where u(a) = 1 for a = 0, . . . , q − 1. (As one can check the corresponding tensor is a fixed point also for the tensor network renormalization flow.) The contraction of four tensors S along the four edges of the square would involve four sums. Due to the (four) delta functions in the sum this, however, reduces to one summation. There will be one delta function δ (q) (i + j − k − l) left, which amounts to the Gauß constraint for the effective tensor T ijkl : This new tensor will in general not be of the factorizing form (5.6) anymore, so generically the decomposition into tensors S will now involve q 2 non-vanishing singular values. That is in case D c < q 2 , the approximation sets in. Nevertheless, the block diagonal form of all tensors and matrices involved can be kept through the following iterations.
At a general iteration step we will work with a tensor T abcd with double indices a = (j a , m a ) where j a = 0, . . . q − 1. As will be explained below the range of m a depends on the SVD in the previous iteration step. These tensors will satisfy the Gauß constraints at all iteration steps. Similarly as before we can define for the even vertices the matrix M (ja,ma)(j b ,m b ),(jc,mc)(j d ,m d ) 1 = T (ja,ma)(j b ,m b )(jc,mc)(j d ,m d ) . Due to the Gauß constraint (5.11) this matrix can be brought into block diagonal form, as the non-vanishing entries must verify j a + j b = j c + j d =: i. We will denote these blocks by M 1 (i). Then the singular value decomposition can be applied on the single blocks M 1 (i) for i = 0, . . . , q − 1, so that This yields of course the same singular values as for the entire matrix. Apart from providing a faster algorithm [74,75], the numerical implementation leads also to more stable results for the following reason. Generically one encounters the case of having singular values with multiplicities higher than one. In this case the singular value decomposition is not unique as one can perform rotations among the basis vectors associated to a given singular value. Here, keeping the block structure explicit prevents a mixing between different blocks induced by these rotations.
To implement the approximation we now have to select the D c largest singular values among the singular values of each block. The number N 1 (i) of singular values selected from the block i determines the range of the second index m i in the double index (i, m i ). That is this number N 1 (i) can take values between zero (no singular value selected) and D c (all selected singular values come from one and the same block). Note that for the first iteration step N 1 (i) = 1 (in case all the u(a) are non-vanishing). We define the S matrices by Note that we have i = j a + j b and i = j c + j d (mod q) for the non-vanishing entries of S 1 and S 2 respectively. For the matrices (M 2 (i)) (jc,mc)(j b ,m b ),(ja,ma)(j d ,m d ) at the odd vertices we proceed similarly, which will result in matrices S . For the nonvanishing entries of these matrices we have j b − j c = i and j d − j a = i (mod q) respectively. The range of the indices m i is now determined by the number N 2 (i) of singular values selected from the block i of the matrix M 2 .
The contraction of the four S-matrices along the square now results in a sum Here the range of m c , m k−j−c is determined by N 1 (c), N 1 (k − j − c) from the previous iteration step whereas m i−c , m j+c is determined by N 2 (i−c), N 2 (j +c) respectively, also from the previous iteration step. Accordingly, the range of m i , m k is determined by N 1 (i), N 1 (k) of the iteration step under consideration (that is by the last SVD) and m j , m l by N 2 (j), N 2 (l) respectively. This altered algorithm has not only the advantage of keeping the Gauß constraint explicit but is also keeping any physical interpretation that might be attached to the representation labels / indices j a . For the gravitational spin foam models these would carry information on lengths and area variables -and the procedure here would correspond to a blocking where the microscopic geometrical variables are basically added to obtain the coarse grained variables. The Gauß constraints arise because of reasons rooted in representation theory: namely that for Abelian groups the tensor product of representations k and k leads to a representation k + k . Indeed the tensors T and S can just be seen as intertwining maps between (tensor products of) representation spaces, and the first index j a of the double index (j a , m a ) encodes the representation carried by the associated edge.
The same arguments based on representation theory apply for the non-Abelian spin net models. Also there the tensors T and S are intertwining maps between representation spaces leading to matrices of block diagonal form. For the abstract spin net models the initial tensor T is basically determined by the choice of projectorP v in (2.9), which restricts the intertwining map from one of maximal rank to one of some smaller rank. Here it will be interesting to study whether any of the projector properties are preserved under coarse graining.
Tensor network methods thus have the potential to give direct insight in the behavior of the projectorsP v under coarse graining, which are the key dynamical entities in spin foam models. Moreover, the kind of coarse graining in the spin net or spin foam representation uses the representation theory underlying these models and furthermore corresponds to a geometrically natural coarse graining.

Equivalence of models
A crucial point in every renormalization method is the question which space of models one is considering, that is in which space the renormalization flow is taking place. Often the coarse graining process leads to models outside this space and usually some approximation method is employed to project the models back into the chosen space. For instance, in the Migdal-Kadanoff approach one considers models with local interactions, that is non-local terms have either to be neglected or replaced by appropriately chosen local interaction terms. In particular for a Z q gauge model with local (single plaquette) interactions one always stays in the form of the Z q gauge model (with single plaquette interactions).
In contrast the space of models in the TNR method is defined by the chosen cutoff D c on the index range of the tensors T . Hence different models, for instance spin net models with different Z q groups, can be considered in the same space. Indeed, for certain initial conditions the initial tensors T abcd might actually define the same tensor network models.
Consider, for instance, the Abelian cutoff models for different (even) q but for the same cutoff parameter K. That is we deal with the initial tensor (5.6) where with k running from −( q 2 − 1) to q 2 (we consider even q). As T abcd = u(a)u(b)u(c)u(d)δ (q) (a+b−c−d) the corresponding (symmetric) matrices M ab,cd 1 and M bc,ad 2 will have zero rows and columns if these include an index a with u(a) = 0. In the SVD these rows/columns can just be neglected as it leads to a vanishing singular value. After neglecting zero rows and columns in the matrices arising from different models, these matrices might however coincide (with appropriate matching of indices) and in this sense define the same TNW model.
Indeed, one can check that in the case of the models (5.15) this occurs for a fixed cutoff K but for varying q as long as q ≥ 4K + 2. In this case the first singular value decomposition for the matrices M 1 , M 2 give for both matrices the same 4K + 1 non-vanishing singular values λ max , λ max − 1, λ max − 1, λ max − 2, λ max − 2, . . . , 1, 1 , where λ max = 2K + 1 . (5. 16) As long as the cutoff in singular values D c is chosen such that D c ≥ 4K+1 the first decomposition (5.2) of the T matrices into S matrices is exact.
In particular this means that Abelian spin net models with q ≥ 4K + 2 and fixed K but different q should go through the same renormalization sequence and hence should also end in the same fixed point. We will see in section 5.4 that this holds almost always in the numerical simulations. There will however be also examples where the fixed points depend on q. In these cases the simulations approach an unstable fixed point and the difference in the simulations appear only for large iteration numbers (more than 100 iterations). Hence, the difference can be explained by numerical instabilities and the fact that our method depends on the parameter q in order to keep the block structure of the T matrices explicit: q defines the number of blocks and hence the kind of possible field redefinitions which underly the method.
Disregarding this point the equivalence between models should also hold if we send q → ∞, that is consider U (1). In general, we see that the TNR method might also be applied to Lie groups (which would lead to infinite dimensional matrices) as long as we consider initial data such that the initial matrices reduce to finite dimensional ones due to either the appearance of zero singular values, or singular values which are sufficiently small. This equivalence between models appears only approximately in the Migdal-Kadanoff method, as there the number of parameters on which the renormalization flow acts is fixed by q, or more generally the size of the group, on which the model is based.

Structure of fixed points
In the Migdal-Kadanoff scheme we considered a renormalization flow within a space described by q (or (q − 1) after normalization) parameters for spin nets based on the group Z q . We are much more flexible with the TNR method, where the number of parameters is determined by the chosen cutoff D c on the number of singular values. Hence there are potentially many more fixed points. Indeed the fixed point structure is considerably more complicated, in particular for D c >> q, as we will exemplify below with the Ising model, q = 2. Note that (with the exception of D c = 2, which can be treated analytically) we will only discuss fixed points which we found as a result of the renormalization process, i.e. by flowing to these fixed points. That is we will miss most of the unstable fixed points, which are those with repellent directions and would require a fine tuning of parameters to flow into.
The main feature of the TNR method -with the approximation based on the singular value decomposition -is the appearance of non-isolated fixed points. These are argued [62,104] to be due to short scale degrees of freedom which are not averaged out by the approximation method employed here. In [62] different forms of additional approximation steps (termed entanglement filtering) are suggested, that apply once the flow reaches the non-isolated fixed points. We will not consider these additional steps here and just describe the type of fixed point encountered in the original TNR method. For future work one should, however, study this issue in more detail. In particular, one could benefit from a detailed comparison between approximations of Migdal-Kadanoff type and approximations based on the TNR method and furthermore from an analysis of relevant and irrelevant directions for the linearized renormalization flow around this kind of fixed points.
In the following we will describe the results of the TNR method for the Ising model in the spin net representation (associated to the high temperature expansion) for the choice of different cutoffs D c , see also [62] for a treatment of the Ising model with the TNR method improved by entanglement filtering.
The smallest cutoff which would also allow a representation of the low temperature fixed point (LTF) u(0) = u(1) = 1 is D c = 2. We start with configurations u(0) = 1, u(1) = x. For x > 0.60858 these flow to the LTF represented by u(0) = u(1) = 1. Accordingly, the singular values λ i for this fixed point are λ 0 = 1, λ = 1 with every block j = 0, 1 contributing one singular value. (We normalize the tensors in every step such that the largest singular value is equal to one.) For 0 ≤ x < 0.60857 the configurations flow towards the high temperature fixed point (HTF) represented by u(0) = 1, u(1) = 0. This corresponds to having only one non-vanishing singular value λ 0 = 1. The transition point at u(1) = 0.6085 corresponds to a phase transition temperature of kT c ≈ 2.572. This is a better approximation to the exact result of kT c ≈ 2.269 than the isotropic Migdal-Kadanoff result of kT c = 3.282 in section 4.1.
In this case, D c = 2, one can analytically compute the flow in some suitable parameters, for instance the tensor components T abcd . This allows to find a further (unstable) fixed point, given by the matrix M ab,cd 1 = T abcd (the numbering of the rows and columns is ab = 00, 11, 01, 10) Note that the tensor describing this fixed point is no longer of the original form (5.6) as this would require s 2 = s 2 1 . Nevertheless it leads to the phase transition between the low and high temperature regime.
Let us now turn to the case D c = 4. In this case configurations with x > 0.6466 will still flow to the LTF characterized by two non-vanishing singular values λ 0 = λ 1 = 1, one from every block. For configurations with x < 0.6465 we encounter however a new type of non-isolated fixed points, so-called Corner Double Line tensors [62,104]. Figure 12: Tensor with Corner Double Line structure 'Double Line' indicates that we have to deal with pairs A = (a, a ) of indices. The tensors are defined by associating a matrix C ab to each corner of the four-valent vertex, see figure 12, Hence if the SVD of C ab includes n non-vanishing singular values η I , we will obtain n 2 nonvanishing singular values of the form η I η I for the matrices obtained from T . The fixed point we encounter in the high temperature region for D c = 4 is of CDL type with corner matrix That is in the SVD for the matrices associated to the tensor T we encounter four non-vanishing singular values λ i = 1, y, y, y 2 . Here y ranges from zero and seems to get arbitrarily close to 1 for u(1) = x reaching the transition point at 0.6496562.... Note that here we encounter a continuous fixed point family of 'high temperature type' ranging from having only one nonvanishing singular value equal to one to having four non-vanishing singular values, with three of these almost equal to one. If we take the transition between this fixed point family and the LTF (with two non-vanishing singular values equal to one) as the phase transition we obtain a critial temperature of kT c = 2.221. This again is a better approximation to the exact phase transition temperature than the D c = 2 result. The appearance of these non-isolated fixed points is interpreted [62,104] to be due to short range entanglement, which is not filtered out by the renormalization flow. Hence certain microscopic details of the models are remembered and the usual universality property of phase transitions does not apply. One can nevertheless argue that the correspondence between phases and fixed points does hold: here models in a certain phase flow to a certain type of fixed points. Indeed we will see that for D c = 9 two different types of non-isolated fixed points appear, representing the low and high temperature regime respectively. In [62] additional approximation steps are suggested, designed to filter out this short range entanglement. We will leave the investigation of these and other additional approximation steps for future work.
The appearance of the CDL type fixed points suggest that cutoffs D c = n 2 with n a natural number, might lead to a fast convergence, as these accommodate the n 2 non-vanishing singular values of a CDL type fixed points. We therefore will also discuss the case D c = 9. Here, both the low temperature phase and the high temperature phase are described by non-isolated fixed points. For the initial value u(1) = x ≤ 0.64293 we flow to fixed points, described by nine non-vanishing singular values of the form λ i = 1, y, y, y 2 , z, z, yz, yz, z 2 . where T CDL is a tensor based on the corner matrix (5.21) and T LT F is the tensor associated to the low temperature fixed point, that is of the form (5.6) with u(0) = u(1) = 1. Therefore eight non-vanishing singular values appear, as products of the four singular values λ = 1, y, y, y 2 for the CDL tensor and the two singular values λ = 1, 1 for the LTF tensor. The appearance of this type of fixed point tensors of product form was conjectured in [62]. We will encounter more fixed points of this form in the next subsection. The transition between the high temperature family of fixed points and low temperature family of fixed points corresponds to a critical temperature of kT c ≈ 2.274, which again approximates the exact result kT c ≈ 2.269 better than the D c = 4 result (kT c ≈ 2.221).

Analysis of Abelian cutoff models
We will now discuss the renormalization behavior of Abelian cutoff models to compare it to the results obtained with the Migdal-Kadanoff method in section 4.2. We will consider the same kind of configurations as in section 4.2, parametrized by the size of the group q and by the cutoff 11 K. Now, as was pointed out in section 5.2 configurations with the same cutoff parameter K but different q might encode the same tensor network model. To utilize this we will choose the same cutoff D c (D c = 16, 25 and equal to 32 for some models) for different q, so that the renormalization flow should also be the same in these cases. This is different from the MK method, where the number of parameters always depends on q.  We have seen in section 4.2 that the MK renormalization flow for the 2D spin net models was much more involved than for the 3D spin foam models: For sufficiently large groups Z q configurations would undergo a behavior determined by unstable fixed points and rather weak phase transitions leading to very long convergence times. Therefore we have to expect similar properties to appear in the tensor network renormalization. Furthermore, as is usually the case if configurations approach a region around phase transitions, the approximation implemented by the cuttoff D c might get less and less accurate. Indeed, the values of the neglected singular values might become comparably large (of order 10 −1 of the largest singular value, which will always be normalized to one), even for the quite large cutoff D c = 32. The TNR method, however, offers the possibility to increase the cutoff and thereby study the influence of the cutoff on the results obtained.
Another issue that appears during a number of simulations are configurations for which the symmetry k → −k (equivalent to reversing all the edges) is broken. The reason for this is that singular values usually appear with a two-fold degeneracy, namely one from a k-block, the other from a (−k)-block. The singular values from the block with k = 0 and k = q/2 are an exception to this rule. Also, quite often singular values appear with even higher degeneracy. Now, taking only a fixed number of singular values into account, one will quite often dismiss one singular value of a degeneracy pair and in this way break the edge reversing symmetry. As the number of k = 0 singular values, which are taken into account by the cutoff is not known a priori, and might even change during the iterations, it might be quite difficult to find a cutoff where this issue does not appear.
In regions 'near' phase transitions, that is where convergence takes very long, these nonsymmetric modes typically lead to unstable, oscillating behavior. This also signifies that the influence of the cutoff is not negligible. In the Migdal-Kadanoff method non-symmetric configurations may only appear due to numerical inaccuracies and it is quite straightforward to implement a symmetrization after each iteration step. A similar procedure for the TNR method (which is however not as straightforward to implement as for the MK method) would probably be very helpful to obtain more reliable results as well. Here we will interpret an oscillating behavior over large iteration times as indicating the presence of an unstable or quasi fixed point.
We can broadly summarize the behavior we encountered in the simulations of the Abelian cutoff models into the following classes, wich we also illustrate in figure 13: • The models flow quite fast (typically during 10 to 30 iterations) to the low temperature fixed point or to a fixed point of type T LT F ⊗ T CDL as described in section 5.3. This happens for cutoff models which are not independent of q. The low temperature fixed points come with q non-vanishing singular values equal to unity. The more complicated fixed points with CDL structure appear, in particular, for the higher cutoff D c = 25.
The singular values associated to the CDL structures are however quite small (∼ 10 −5 to 0.2), so that we can safely interpret these fixed points as of low temperature type. It also happens that initially the configurations flow to a T LT F ⊗ T CDL fixed point with a higher number of non-vanishing singular values than provided by the chosen cutoff. In this case the configurations slowly keep changing such that the singular values associated to the CDL structure decrease. With the exception of two models (K = 3, q = 12 and K = 4, q = 16) such a behavior appearing in the D c = 16 simulation would be confirmed in the D c = 25 treatment.
• An interesting special case is K = 1 and q ≥ 6 (which describes the same tensor network model independent of q as long as q ≥ 6). For D c = 16 the fixed point (reached after about 60 iterations) is of the form T LT F ⊗ T CDL , but the T LT F factor comes with only four singular values equal to unity. The CDL structure leads to 16 non-vanishing singular values, where apart from the four singular values equal to one, the next eight ones are equal to 0.86 and the final four equal to 0.73.
The D c = 25 results are slightly different: due to the appearance of non-symmetric configurations as described above, the configurations are oscillating for more than 80 iterations. Then the results do actually start to depend on q: for q = 6 a T LT F ⊗ T CDL fixed point is reached with six singular values equal to one (so it is really a q = 6 low temperature fixed point), 12 singular values equal to 0.81 and six equal to 0.655. Whereas for q = 8 a fixed point is approached (around 200 iterations) with eight singular values equal to one, another sixteen are around 0.3, plus eight singular values around 0.12. This gives more than 25 non-vanishing singular values, so as described above the configuration is slowly changing, decreasing the CDL singular values. Note however that due to the high iteration numbers involved and the appearance of 'non-symmetric' configurations the results should be taken with some care.
• For the D c = 16 simulations there are two examples K = 2 and q = 10 or higher, and K = 3 and q = 14 or higher, which show stable behavior and seem to approach a (nontrivial) fixed point for long iteration times (approx. 20 iterations for K = 2 and approx. Note that examples flowing to a high temperature fixed point do not appear. The best candidate would be the K = 1 configurations (describing the same tensor network models starting with q = 6 and larger q), which however flow for D c = 16 to a fixed point with four singular values equal to unity, plus more non-vanishing singular values due to the CDL structure.
In the following  (14) QF (6)→Osc (13) QF (8)→Osc (53)  2 12 LT (59) QF (6)→Osc (13) QF (8) Table 1: Summary of the renormalization flow of Abelian cutoff models. To a precision of 10 −10 , the numbers in brackets give the iteration times (halfed for TNR to make comparable) it takes to reach a certain behaviour: LT (low temperature fixed point, bold numbers denote the number of ones that appear explicit), HT (high temperature fixed point), QF (quasi fixed point), Osc (oscillating behaviour), CDL (corner double line), NULL (finite part that dies off). For a given K and growing q, the fine line indicates the beginning of physically equivalent models in the sense of section 5.2.
The overall picture of the Migdal Kadanoff results is confirmed by the tensor network simulations: Most configurations flow to the low temperature fixed point or to a low temperature fixed point embellished with a CDL structure. With growing q, cutoff models with sufficiently small K show long stable phases where a fixed point seems to be approached but then enter an unstable phase with (slightly) oscillating behaviour. For some cases these converge finally to a LTF fixed point, however this result has to be taken with some care, due to the oscillating phase in which 'non-symmetric' configurations appear.
There is one important difference between the MK and the TNR results. As we explained in section 5.2 examples with q ≥ 4K + 2 should encode the same physical models both exactly and in the approximation provided by choosing the cutoff D c . This property is inherent in the TNR method (and only appears to be violated for large iteration numbers due to the oscillating behavior). But in the MK method the truncation is determined by q and hence the renormalization flow is still different for q not much larger than q K := 4K + 2. For much larger q the renormalization sequences will turn out to be almost the same, however this behavior sets in later than in the TNR examples. Moreover from [92] we have to expect that going to q → ∞ the configurations should flow to the high temperature fixed point. This is the reason why the quasi fixed point behavior has to set in earlier for the TNR method, as here the results starting with q K should in principle also hold in the limit q → ∞. From this point of view it is interesting that we have not seen an example in the TNR method (especially not K = 1, as opposed to the MK K = 1, q = 8 example) which would flow to the high temperature fixed point. Hence the TNR method might be able to detect two different phases for the U (1) theory.
As we have also seen definite conclusions are very much hindered by the appearance of nonsymmetric configurations, i.e. where the k-blocks would differ from the −k-blocks. These also appeared in the MK simulations for the examples which go through very long almost stable phases. In this case non-symmetric configurations only appeared due to numerical errors and this problem can be easily cured in the MK method by symmetrizing the Q(k) parameters after each coarse graining step.
In the TNR method non-symmetric configurations not only appear due to numerical errors but are also caused by the cutoff, which might neglect one of a pair of degenerate singular values. This happens quite generically if the cutoff is increased (apart from the D c = 25 simulations we also tried D c = 32 and inbetween values) and it is increasingly difficult to find a value for D c where this does not appear, say, for the first twenty iterations.
This non-symmetric behavior does not matter so much for configurations which would converge fast to some stable fixed point but do cause long sequences of oscillating behavior for configurations, which we suspect would otherwise rather approach slowly an unstable fixed point. For these examples the ordered sequence of singular values at a given iteration decreases rather slowly. Even for D c = 32 the neglected singular values can be of order 10 −1 and an unsymmetric cutoff will have considerable influence on the overall behaviour. This problem is especially pronounced for the 2D models with larger q, as we there encounter rather weak phase transitions.
For future work one should address this issue 12 . There are different possibilities, one is to use a symmetric parametrization in the coarse graining procedure, i.e. only work with the k-blocks where k ≤ q 2 . Another option is to use an adaptive cutoff D c for each iteration step, such that it avoids to cut between (k, −k) pairs.
A third option would be to change the approximation scheme slightly by choosing a block dependent cutoff D c (k). This would actually simplify the algorithm considerably. For sufficiently high D c and D c (k) these schemes should yield equivalent results, this has however to be tested.
Note that non-symmetric configurations can however also appear due to numerical instabilities, so that one might have to implement some symmetrization procedure (if one does not work with a parametrization that only allows symmetric configurations).

Discussion and outlook
Extracting large scale physics from spin foam models is one of the most pressing issues for the field. Here we advocated the development and use of coarse graining and renormalization techniques so that hopefully a feasible method for 4D models can be found. To this end, we introduced a wide range of simplified models, spin foams with finite groups and moreover spin nets, which can be seen as a dimensionally reduced version of spin foams. This not only enables us to test and develop coarse graining methods but also to obtain some physical insights and to put forward conjectures about the dynamical behavior of the full models.
This strategy allows us to adopt coarse graining methods from lattice gauge theory and condensed matter system, here the Migdal-Kadanoff scheme and the tensor network renormalization method. The two methods have different drawbacks and advantages. The Migdal-Kadanoff scheme facilitates quick results (even a large number of iterations can take only seconds on a PC), so that an overview of the phase structure encoded in the model can easily be obtained.
Here the main question is how this method can be generalized to non-Abelian spin foams with non-trivial projectors, which no longer fall into the class of standard lattice gauge theories.
The tensor network renormalization method has the advantage of providing a systematic improvement on the accuracy of the results. The required effort is considerably larger (100 D c = 32 iterations may take several days on a PC). The method is however very general and in its version based on the spin net representation, allows direct access to the behavior of the (vertex) projector under coarse graining. Moreover, the blocking of variables is very natural if one takes into account the geometric meaning of the representation labels in the gravitational spin foam models. We presented an algorithm in which the Gauß constraint are kept explicitly intact.
The tensor network renormalization method is easily generalizable to models with non-Abelian groups, and indeed would nicely interact with the group theoretic content of these models, see also [74]. Here it will be very interesting to study how a non-trivial vertex or edge projector might change the phase structure as compared to the standard choice of the Haar projector. This will also facilitate a better understanding of the dynamics in the full gravitational models. To this end, a class of finite group models emulating the current gravitational EPRL models [8] is constructed in [106]. An important question for future research will be whether the degenerate phase, that is the high temperature fixed point can be avoided by selecting suitable projectors. This can already be studied for the 2D spin net models. Furthermore it has to be explored how the tensor network renormalization method can successfully be applied to three and four-dimensional spin net and spin foam models.
An alternative to using the tensor network formalism to describe the partition function of a system would be to apply tensor network renormalization as a kind of improved mean field approach in a canonical quantization. This would change the (statistical) systems from Ddimensional classical to (D − 1)-dimensional quantum ones. In this case the tensor networks would provide an ansatz for the ground states of the models [73,76,105]. In gravity, instead of minimizing one Hamiltonian or energy functional, one rather has to deal with a number of constraints, which have to annihilate the so-called physical states. The master constraint [107][108][109] or the uniform discretization approach [110,111] provide a framework where just one (master) constraint has to be minimized.
In this work we have considered systems on a regular lattice, as this made explicit simulations feasible. Nevertheless one should consider generalizations of these methods to random lattices [49,50]. Another question is how the phase structure found on a fixed lattice relates to phases in models, that include a sum over all possible lattices [112][113][114][115], in particular regarding the remarks in section 3.
We see this work as a contribution towards a closer link between the quantum gravity and the statistical physics / condensed matter communities. For quantum gravity researchers, this offers the prospect of new concepts and numerical tools to study the large-scale physics of their models, whereas for people working in statistical physics it offers new models, new questions, a geometrical perspective and a rich set of mathematical tools behind it.