Coarse graining spin foam quantum gravity -- a review

In quantum gravity, we envision renormalization as the key tool for bridging the gap between microscopic models and observable scales. For spin foam quantum gravity, which is defined on a discretisation akin to lattice gauge theories, the goal is to derive an effective theory on a coarser discretisation from the dynamics on the finer one, coarse graining the system in the process and thus relating physics at different scales. In this review I will discuss the motivation for studying renormalization in spin foam quantum gravity, e.g. to restore diffeomorphism symmetry, and explain how to define renormalization in a background independent setting by formulating it in terms of boundary data. I will motivate the importance of the boundary data by studying coarse graining of a concrete example and extending this to the spin foam setting. This will naturally lead me to the methods currently used for renormalizing spin foam quantum gravity, such as tensor network renormalization, and a discussion of recent results. I will conclude with an overview of future prospects and research directions.

In quantum gravity, we envision renormalization as the key tool for bridging the gap between microscopic models and observable scales. For spin foam quantum gravity, which is defined on a discretisation akin to lattice gauge theories, the goal is to derive an effective theory on a coarser discretisation from the dynamics on the finer one, coarse graining the system in the process and thus relating physics at different scales.
In this review I will discuss the motivation for studying renormalization in spin foam quantum gravity, e.g. to restore diffeomorphism symmetry, and explain how to define renormalization in a background independent setting by formulating it in terms of boundary data. I will motivate the importance of the boundary data by studying coarse graining of a concrete example and extending this to the spin foam setting. This will naturally lead me to the methods currently used for renormalizing spin foam quantum gravity, such as tensor network renormalization, and a discussion of recent results. I will conclude with an overview of future prospects and research directions.

I. A BRIEF INTRODUCTION TO SPIN FOAM QUANTUM GRAVITY
Spin foam quantum gravity [1,2] is a promising approach to quantum gravity closely related to loop quantum gravity [3]. The aim is to define the path integral for gravity in a non-perturbative and background independent fashion, that is without any reference to a fixed background space-time or structure.
The starting point of spin foam models is the Plebanski-Holst formulation of general relativity [4], in which gravity is formulated as constrained topological BF theory [5]. To formulate this theory as a path integral, one introduces a lattice as a regulator, more precisely a 2-complex, in order to truncate the number of degrees of freedom. On this 2-complex, which is a collection of vertices, edges and faces, the topological theory is first discretized and quantized. This is in close analogy to 3D (topological) gravity, where this formulation gives rise to the Ponzano-Regge model [6][7][8], a well-defined model of 3D quantum gravity defined on a triangulation.
However, gravity in 4D is not topological, which requires the implementation of so-called simplicity constraints. In the continuum they serve the role to break the too many symmetries of the theory and reduce the B-field in BF theory to a simple 2-form, reducing the action to the familiar Holst action [9]. In spin foam quantum gravity, one derives such constraints for the discretisation of the classical B-field, so-called bivectors. In 4D, bivectors are assigned to 2D faces, e.g. triangles, and encode their geometry. The constraints ensure that these bivectors are simple, i.e. they can always be written as a wedge product of two vectors. Geometrically these vectors span two edges of a triangle. Different versions of these discrete constraints agree for single, classical building blocks, e.g. a 4-simplex, such that they correspond to different discretisations. However, their implementation in the quantum theory, which leads to restrictions on the variables of the theory, generically results in different models with starkly different dynamics. Two examples are the Barrett-Crane (BC) model [10,11] and the Engle-Pereira-Rovelli-Livine / Freidel-Krasnov (EPRL / FK) model [12][13][14]. The former strongly implements a condition on bivectors, which significantly reduces the degrees of freedom of the model. This was criticized [15,16] and motivated the development of the EPRL / FK model, in which constraints are implemented weakly, i.e. at the level of expectation values. Despite these insights, in remains an open question whether these constraints are sufficient to recover general relativity in a continuum limit. The hope is that coarse graining / renormalization can shed a light on this intriguing question.
Despite these differences, all spin foam models are written in a similar form. The 2-complex, which is frequently dual to a triangulation, is coloured with group theoretic data: each face f carries an irreducible representation ρ f of the underlying symmetry group (Spin(4) for Riemannian, SL(2, C) for Lorentzian signature), while each edge e carries an intertwiner ι e , an invariant tensor in the tensor product of representation spaces associated to faces meeting at an edge. These data encode the geometry of the spin foam: in 4D, each edge is dual to a 3D polyhedron, which has as many faces as (dual) faces that share this edge. Then, the areas of these faces are given by the associated representations ρ f . However, this does not determine the shape of the polyhedron uniquely, see e.g. a tetrahedron. A flat tetrahedron is uniquely determined by specifying its six edge lengths, whereas it only has four faces. Thus, the areas of these faces alone do not fix the shape of the tetrahedron. Part of the information on the shape is stored in the intertwiner, which can be expanded into an orthonormal basis using group representation theory. For a tetrahedron, its dual 4-valent intertwiner can be split into two 3-valent ones, where the new link carries again a group representation labeling the basis element. Geometrically this representation gives the area of a parallelogram spanned by the midpoints of the edges of the tetrahedron [5,17], see also fig. 1. However, due to the uncertainty principle the shape cannot be fully specified since area operators associated with intersecting faces do not commute and thus cannot be diagonalized simultaneously. Note that coherent intertwiners can be defined that are sharply peaked on the geometry of a classical polyhedron 1 . From these 3D polyhedra a 4D geometry is built at the vertices of a 2-complex. At such a vertex, several edges and faces meet, indicating how 3D polyhedra are glued together to form a 4D geometry. If two edges meet at the same vertex and are part of the same face, their dual 3D polyhedra are glued along the shared face. Crucially, since the representation associated to the face determines its area, it is ensured that the face has the same area in both polyhedra. From the group theoretic data, this 'gluing' is performed by contracting the intertwiners according to the combinatorics of the 2-complex, which essentially amounts to a spin network evaluation 2 . The resulting number is known as the vertex amplitude A v , i.e. the amplitude of the spin foam model assigned to the discrete 4D geometry dual to the vertex v with configuration {ρ f , ι e }. Similarly, we assign local amplitudes to the A e and A f to the edges e and faces f respectively. The former ensures that intertwiners are normalized, while the latter corresponds to the dimension of the representation ρ f . See [1] for more details of the derivation. Eventually, the path integral is defined as a sum over all these configurations: Crucially, these geometric building blocks and amplitudes are derived from general relativity formulated as a constrained topological field theory. In case the spin foam has a boundary (see fig. 2), it serves as an amplitude functional mapping states from its boundary Hilbert space into the complex numbers. This concept will be crucial in this review. At the level of a a few simplices, these models are well explored. A well understood result across models, which furthermore underlines the relation to general relativity, is the asymptotic expansion of the vertex amplitude dual to a 4-simplex [21][22][23][24][25]. In these works the vertex amplitude is investigated for coherent intertwiners, which are sharply peaked on the geometry of classical polyhedra. Then the vertex amplitude can be written as an integral over several copies of the symmetry group. This integral is then evaluated via a stationary phase approximation by uniformly scaling up all representations. Hence it is commonly referred to as large-j limit, in which the amplitude is generically proportional to the cosine of the Regge action [26,27], a discretisation of classical gravity. Physically, this amplitude should be valid for 4-simplices of mesoscopic and even macroscopic size. In recent years, numerical calculations of the vertex amplitude beyond the asymptotic expansion have seen promising progress [28,29].
FIG. 2. A spin foam in 3D with boundary. The boundary of this spin foam is made up of an initial and a final graph, which carry states ψi and ψ f respectively. The links of a graph carry representations ρ f , whereas the nodes carry intertwiners ιe. For ψi, these are ρi, i = 1, 2, 3, and ι1. The data for ψ f is omitted for clarity of the figure. The 2-complex interpolates between these two graphs and evolves the states. Thus, in the 2-complex, representations are associated with faces and intertwiners with edges. Several edges meet at a vertex v, here shown in blue. The vertex shown here is dual to a 3D tetrahedron: four edges, each dual to a triangle, are glued together to form a tetrahedron.
Beyond a single building block, the dynamics is less explored, in particular how the choice of the 2-complex impacts the results of the theory. Indeed, a priori the theory itself does not specify how the 2-complex should be chosen. In this review we take a viewpoint that is akin to lattice gauge theory, and regard the 2-complex as a regulator, a particular choice to truncate the number of degrees of freedom of the theory. As such, physics must not depend on this choice and it is must be removed eventually, e.g. in a continuum limit, in order to derive consistent results. One route towards such a limit lies in coarse graining: By coarse graining, i.e. defining an effective coarse amplitude from a collection of fine ones, we readily relate two theories defined on two different regulators. Moreover, by coarse graining we gain insight into the dynamics of a collection of building blocks and learn which configurations are more relevant on a coarser scale. In short, the aim is to derive a family of amplitudes to assign to different regulators, which reproduce the same physics (at least approximately). This defines a renormalization group flow of amplitudes [30,31]. It is the purpose of this article to review the progress of this approach and outline how it helps turning spin foam quantum gravity into a computational formalism.
This review is structured as follows: in section II we start by outlining the most pressing challenges faced by spin foam quantum gravity and how these are addressed by coarse graining. Section III discusses the issue of restoring diffeomorphism invariance in the discrete as well as the typical appearance of non-local interactions under coarse graining, which is one motivation for the consistent boundary formulation outlined in section IV. Section V reviews two numerical methods to perform such coarse graining algorithms, namely tensor network renormalization and restricted path integral models. In section VI we conclude with several interesting future research directions.

II. KEY CHALLENGES IN SPIN FOAM QUANTUM GRAVITY
Before explaining renormalization in spin foam models and its progress over the last decade, it is crucial to first discuss the key challenges spin foam quantum gravity is facing and how renormalization plays a vital role in overcoming them.

A. Fate of diffeomorphism symmetry
Diffeomorphism symmetry, as the fundamental symmetry of general relativity, is deeply intertwined with the dynamics of gravity. It implies that physics must not depend on the choice of coordinates and only diffeomorphism invariant (Dirac) observables are physically meaningful [32,33]. Moreover, this symmetry forbids a choice of a preferred or fixed background space-time. Conversely, the complexity of this symmetry is a root of the difficulty for defining a theory of quantum gravity; spin foam models are no exception.
While spin foam quantum gravity embraces the concept of background independence, the introduced regulator, frequently a 2-complex dual to a triangulation, generically breaks (a discrete remnant of) diffeomorphism symmetry [34][35][36], often called a vertex translation symmetry [37]. There exist instances, where this symmetry is preserved in the discrete, where the discretisation perfectly reflects the continuum dynamics, or the symmetry can be restored iteratively via coarse graining. We explain this in detail in section III. For spin foams to be a viable theory of quantum gravity, diffeomorphism symmetry must be restored, at least approximately, in order to derive reliable physical predictions. There exists strong evidence that the amplitudes of the system can be systematically improved via coarse graining [38][39][40], such that the symmetry is broken less. Naturally the question arises whether this procedure converges to a fixed point, which would automatically imply an independence of the chosen regulator. Due to the non-local nature of diffeomorphism symmetry and in order to find a theory with propagating degrees of freedom, we conjecture such a fixed point to lie on a phase transition of second order. There it would be possible to take the continuum (or rather refinement limit) of the theory.

B. Discretisation (in)dependence
Closely related to diffeomorphism symmetry is the question of discretisation (in)dependence. Generically the results computed in spin foam models will depend sensitively on the chosen regulator, e.g. the number of simplices and subsequently the number of degrees of freedom. Moreover, there is no input from the theory itself which regulator to choose. However, in order to have a viable theory, it is imperative to find the same results no matter which discretisation is chosen, at least to an approximation.
In the research community, there exist two complementary paths addressing this question [41]. On the one hand, there is the approach to solve discretisation dependence by summing over all possible regulators, e.g. triangulations. This summation over triangulations (and topologies) is most holistically formulated in terms of group field theories [42,43], which are quantum field theories formulated on several copies of a Lie group. The fields represent atoms of space-time, e.g. tetrahedra, whose interaction terms describe how 4D objects are formed, e.g. five tetrahedra glued together to form a 4-simplex. From this formalism, spin foam amplitudes arise as Feynman diagrams in a perturbative expansion. As for all quantum field theories, it must be shown that this theory is renormalizable, e.g. via perturbative or non-perturbative methods, see [44] for a review.
On the other hand, we discuss the refinement approach [30,31] in this review, where we interpret the triangulation as a regulator to truncate the number of degrees of freedom, similar to the lattice in lattice field theories. The idea to overcome discretisation dependence is by assigning different amplitudes to different discretisations in such a way that the results agree. One example is to derive coarse amplitudes from fine ones via coarse graining. In this way, we are relating theories across different discretisations. The goal is to derive such relations for all possible discretisations, which is equivalent to a complete renormalization group trajectory. Again, this is similar to lattice field theory, where one also assigns different theories to different lattices, parametrized by coupling constants.

C. Computability
The choice of a discretisation (and thus appropriate) amplitude also enters, at least partially, in another key challenge for spin foams, namely their computability. To be more precise, by computability we refer to two interconnected issues. On the one hand, there is the challenge to compute the fundamental spin foam amplitudes for a single building block, e.g. a 4-simplex. While this is well-studied and explored in the semi-classical regime [21][22][23][24][25], in particular using coherent states and stationary phase approximation, computing a vertex amplitude in the quantum regime, e.g. for small spins, can only be done numerically. However, in recent years there has been significant progress in computing these amplitudes, e.g. for the EPRL / FK in Euclidean and the more challenging Lorentzian signature [28,29,45].
Renormalization and coarse graining become important at the stages when we calculate amplitudes or observables for mulitple vertices / larger triangulations. Even if we have an efficient way of calculating spin foam amplitudes (or can access the relevant amplitudes from a database), summing over the various degrees of freedom remains a difficult task for such a high dimensional configuration space 3 . However, if we assume that the full RG trajectory of the system is known, we can use the discretisation to our advantage and perform the same calculation on a much coarser spin foam with appropriately adapted amplitudes. Alternatively and more realistically, one can envision coarse graining the system first, essentially evaluating it in parts, deriving an effective theory on a coarser regulator from a finer one. On this coarse theory, expectation values of coarse observables can efficiently computed. This method is already realized nowadays in tensor network renormalization techniques [46][47][48], see e.g. [49]. Note that the existence of a continuum limit is not assumed, rather we assume that coarse graining can be performed without severe truncations.

D. Uniqueness, phase diagram and continuum limit
Discretizing a continuum theory is generically not a unique process, take the 1D non-relativistic particle in a nontrivial potential as an example. There exist many choices how to discretize the potential, which all result in different dynamics. However, the expectation is that, no matter the choice, we reobtain the original continuum physics in a suitable continuum limit (or approximate it well in a fine discretisation). This is even more severe in the case where the continuum theory possesses a symmetry, like reparametrisation or diffeomorphism invariance, which is broken in general in the discrete [37,38].
These topics, uniqueness of the theory, universality and the continuum limit remain open questions in spin foam quantum gravity. Modern spin foam models are frequently derived by starting from topological BF theory and then imposing simplicity constraints in the discrete [1]. The latter procedure is not unique, where e.g. the well-developed EPRL / FK model imposes the linear simplicity constraints weakly [13,14]. Some effects on different choices of (implementations of) simplicity constraints can be found in the literature [15], however a phase diagram differentiating different universal dynamics is missing, and with it potential hints for a continuum limit and UV-completion of the theory.
These key challenges are deeply intertwined with one another and can be addressed by a coarse graining / renormalization scheme. In the following we review how our understanding of these connections developed over time and what the role of coarse graining is.

III. RESTORING DIFFEOMORPHISM SYMMETRY IN THE DISCRETE
Regge calculus [26] is a discretisation of general relativity. In it the differentiable manifold is replaced by a Ddimensional triangulation, whose edge lengths are the dynamical degrees of freedom. Crucially, Regge calculus does not refer to coordinates of vertices of the triangulation and is solely formulated in terms of their distances. Hence it is manifestly coordinate free. Each of the D-simplices is internally flat, i.e. its D + 1 vertices can be embedded into R D . Curvature is distributional and located on (D − 2)-sublimplices, so-called hinges. To each of these hinges in the bulk one associates a deficit angle h , which is the difference between the sum all dihedral angles of simplices meeting at this hinge and 2π. This is nicely visualized in d = 2: Several triangles meet at a single vertex. If their angles located at this vertex sum up to 2π, it is flat and can be drawn on a piece of paper. However, if the deficit angle differs from 0, e.g. if h > 0, we can no longer embed this collection of triangles into R 2 and observe positive curvature around that vertex. Note that the edge lengths are the only dynamical variables, as the dihedral angles are given as functions of the edge lengths.
In addition to making no reference to coordinates, in some instances Regge calculus possesses additional symmetries in the discrete linked to diffeomorphism invariance [50]. One such example is 3D Regge calculus for Λ = 0: its equations of motion state that all deficit angles e = 0 in the bulk, for all boundary data, describing a theory that glues piecewise flat tetrahedra in a flat way. Thus, it perfectly matches the continuum solution. Moreover, the Regge action is invariant under vertex translations, i.e. moving a vertex and accordingly changing the edge lengths it is connected to. One such example is the 4-1 Pachner move: If we place an additional vertex in the centre of a tetrahedron, we can freely choose three edge lengths connecting it to the vertices of the coarse tetrahedron. The fourth is then fixed uniquely by the equations of motion. This symmetry is reflected by nulleigenvalues of the matrix of second derivates of the action.
Moreover, the 3D Regge action itself is invariant under such Pachner moves, i.e. local changes of the triangulation. This renders it triangulation independent, since any triangulation of a manifold can be related to any other triangulation of the same manifold by a consecutive application of Pachner moves [51,52]. This is not surprising, since 3D gravity is topological, i.e. has no local degrees of freedom. Nevertheless, we are convinced that triangulation independence and diffeomorphism symmetry in the discrete, in the form of a vertex translation symmetry, are closely related also beyond topological theories. Diffeomorphism symmetry is deeply entangled with the dynamics of general relativity. When perfectly realized in a discrete system, by fully capturing the continuum dynamics, it is irrelevant whether we consider a coarse or a fine discretisation. Thus the theory is discretisation independent. Invoking the invariance under vertex displacements, we can imagine this by moving vertices on top of each other, effectively removing them. Conversely, achieving discretisation independence by finding a fixed point of a coarse graining flow, e.g. on a second order phase transition, does not necessarily imply that diffeomorphism symmetry is restored, yet this conjecture is supported by several examples that we outline below.
A nice example how coarse graining can improve an action (or amplitudes in the quantum case) is again 3D Regge calculus with a non-vanishing cosmological constant. Due to the cosmological constant, the equations of motion state that deficit angles do not vanish. Moreover, the theory is not triangulation independent and the vertex translation symmetry is broken. In [38], Bahr and Dittrich device a coarse graining scheme for the triangulation: On a refined triangulation, subdividing large edges into small ones, they solve the equations of motion for the small edges and define an effective action for the remaining large ones. This procedure converges to a fixed point action, which describes Regge calculus for constantly curved tetrahedra. On this fixed point, deficit angles vanish, the theory is triangulation independent (by definition) and the vertex displacement symmetry is restored. Indeed, this improved discretisation / action encodes the continuum solution in the discrete, thus implementing a discrete remnant of diffeomorphism symmetry. Moreover, since it correctly captures the continuum dynamics, no information or accuracy is lost when using coarse triangulations. An analogous quantum version is the Turaev-Viro spin foam model [53], defined as a quantum deformed Ponzano-Regge model [6].
There exist several instances where the continuum solution is pulled back to the discrete setting, where the discrete theory possesses a vertex translation symmetry. One example is 4D Regge calculus [50], when the boundary data allow for flat solutions in the bulk, or the 1D quantum parametrised (an)harmonic oscillator [37]. In general we cannot guess these solutions, but with coarse graining methods we can construct or at least approximate them well. However, the examples that we discuss here are either topological or one-dimensional, and hence it is possible to retain a local description. For higher dimensional, interacting theories, non-local interactions appear, which can be a stumblestone for coarse graining methods.

A. Non-localities
Before explaining non-localities or non-local interactions, we first need to state what a local theory is in this context. In most discrete theories, we associate variables to parts of the discretisation, e.g. in spin foam models we assign irreducible representations ρ f to faces f of the dual complex and intertwiners ι e to edges e. In Regge calculus, we assign edge lengths l e to the edges e of a triangulation. We define this theory to be local if the partition functions is given by a product of amplitudes assigned to (sub)simplices or if the action is given as a sum over actions assigned to (sub)simplices. Moreover, the action and amplitude for each (sub)simplex only depend on those variables attached to (sub)simplices contained in the (sub)simplex. Spin foam models are an example for such local theories, since the partition function is given via a local assignment of vertex, edge and face amplitudes, see eq. (1). Similarly, the Regge action can be written as a sum over contributions associated to the D-simplices of the triangulation.
When we apply coarse graining methods to such interacting, i.e. non-topological, theories, it is highly unlikely that this local form of the theory can be preserved. There exist several examples in the literature where this has been shown in the past. In [54], 4D Regge calculus was linearized around a flat background solution and the perturbations of the edge length integrated over. The question is whether it is possible to find a path integral measure that is invariant under Pachner moves. However, when integrating out these degrees of freedom, one picks up a non-local factor that cannot be written as a local product. In [55] it is shown that said factor is related to a condition whether the six vertices involved in the Pachner move lie on a 3-sphere. Moreover, these articles reveal that the 4D Regge action itself is not invariant under Pachner moves. In a similar vain, [56] studies Pachner moves in 4D holomorphic spin foam models [57]. The advantage of these models is that Pachner moves can be explicitly computed. Again, the resulting amplitude is non-local, in the sense that the resulting expression cannot be written as a assignment of local amplitudes.
To illustrate this point further let us consider the concrete example of the 2D Ising model.

Ising model as an example
There exist plenty of ways to coarse grain discrete systems. A straightforward example is the 2D Ising model subject to a simple decimation procedure, where one simply sums over "every other" spin to derive an effective model on a larger scale.
We consider the Ising model defined on a regular 2D lattice with vanishing external magnetic field. There are only nearest neighbour interactions, i.e. an Ising spin σ i ∈ {−1, 1} only interacts with its direct neighbours. Then we can write the partition function as product of weights associated to the edges of the lattice: where β is the inverse temperature, and s(e) / t(e) denote the source and target of the edge e 4 . Note that the system has a global Z 2 symmetry, it remains invariant if all spins are flipped.
We implement a decimation procedure by summing over every other spin, essentially evaluating the partition function in parts. In order to derive the new effective amplitude of the system, it is sufficient to consider four Ising spins σ 1 , . . . , σ 4 that all connect to another Ising spinσ, see fig. 4. The four coarse spins sit on the corner of a coarser square rotated by 45 • with spinσ in the center of the square. We obtain: σ exp(βσ (σ 1 + σ 2 + σ 3 + σ 4 )) = 2 cosh(β (σ 1 + σ 2 + σ 3 + σ 4 ). ( Clearly this expression is not of the same form as the original action, in particular is not written in terms of Z 2 group multiplications. Nevertheless the remaining spins still satisfy the global Z 2 symmetry. Thanks to this global symmetry, this expression can only take three different values depending on the configurations of the four spins {σ i }; either all spins are aligned, one is not aligned with the others or we have two pairs of aligned spins. To express this again in terms of spin interactions, we make the most general ansatz of four spin interactions compatible with the global Z 2 symmetry: a is the parameter for nearest neighbour interactions, b for next-to-nearest neighbour interactions, c for a four spin interaction and d is a constant. We can compare these equations directly for each configuration: where we denote a cyclic order i, j, k, l around the square. Here we have four equations for four unknown parameters, which we can straightforwardly solve. We leave deriving the solution to the interested reader. The coarse grained amplitude is notably different than the initial one. While we find again nearest neighbour interactions, new non-local interactions appear as well. From this new form it is not obvious how to return to the original expression. Moreover, it is not clear how to iterate the procedure without approximations: decimating one spin alone results in non-local interactions among eight spins, some of which ought to be decimated as well, see again fig. 4. Nevertheless, already this simple example hints towards a resolution: the non-localities arise since we attempt to express the coarse grained dynamics in terms of the old degrees of freedom and building blocks. Yet we can still write the partition function as a local product of amplitudes associated to rotated squares, where the non-local interactions are completely contained within these locally assigned amplitudes. In the next section we introduce this change of perspective more concretely and discuss the concept of generalized boundary data 5 .

IV. CHANGE OF CONCEPT: GENERALIZED BOUNDARIES AND AMPLITUDE MAPS
The vital insight to arrive at a practical coarse graining scheme for spin foams is the following: instead of pertaining the original degrees of freedom and building blocks, e.g. simplices, and allowing more and more complicated, non-local interactions among them, we work with locally interacting amplitudes, which allow for more general and complex boundary data. The non-locality is still present, yet contained within the amplitudes and expressed as interactions of these boundary data. Thus, the complexity of the boundary data controls the non-locality preserved under coarse graining and the complexity of the amplitude. Truncating the boundary data allows us to introduce controllable approximations, while the partition function is still written as a local assignment of amplitudes. Thus we can iterate a coarse graining procedure that only needs to consider few building blocks at a time.
As a path integral approach, spin foam quantum gravity is already phrased in this language, as amplitude functionals for certain boundary states. Take a spin foam on a 2-complex Γ b with boundary b. Since the 2-complex is discrete, namely a collection of faces, edges and vertices, its boundary b is also discrete, namely a graph, with nodes and links. The complexity of this boundary depends on the number of nodes and links. To each of these boundaries b one associates a boundary Hilbert space H b , whose complexity again depends on the complexity of the boundary. A spin foam model for said two complex then acts as an amplitude functional A b mapping states The vital difference is that we allow for more general building blocks, in particular with more complex boundary data and thus boundary Hilbert spaces. When using Pachner moves, one integrates out bulk degrees of freedom while keeping the boundary unchanged. In 4D, when performing a 4 − 2 move, one integrates out bulk variables and derives one effective amplitude for two glued 4-simplices, prescribed by the same boundary data. However splitting this effective amplitude into two, one assigned to each building block, is not straightforward due to the previously mentioned non-local interactions. Instead, we allow for more general building blocks with more complicated boundary data. That way, we still have local assignments of amplitudes to building blocks, which in turn interact locally with neighbouring ones. In turn, non-localities still arise, yet they are contained in each building block and captured by more complex boundary data.
While this picture recasts the problem of arising non-localities, three immediate challenges arise. Firstly, iterating this procedure leads to more and more complicated building blocks, whose amplitudes are still given by the fine boundary degrees of freedom. From this perspective we have not achieved a derivation of coarse scale physics, since the dynamics are still expressed in terms of fine scale degrees of freedom. Secondly, in order to define a renormalization group flow it is crucial to compare amplitudes after each coarse graining step. And thirdly, deriving amplitudes with more and more complex boundary data quickly becomes unfeasible, independent whether one is using analytical or numerical techniques, as one can already see for the 2D Ising model.
Hence, the next vital ingredient for a coarse graining scheme is the introduction of variable transformations, that map a collection of fine boundary degrees of freedom to a collection of coarse effective degrees of freedom. More precisely, we want to map states on a fine boundary Hilbert space H b on b to a coarse boundary Hilbert space H b on b. In the next section, we will explain the idea behind this concept and its physical interpretation. To do so, we work in the opposite direction and explain how to add degrees of freedom using embedding maps.
A. Embedding maps and the notion of vacuum As outlined above, a key idea of any renormalization procedure is to compare and relate theories defined on different scales. Given two spin foam amplitudes A b and A b , which are functionals for the Hilbert spaces H b and H b respectively, these amplitudes can only be compared for the same physical processes. That is, given a state ψ b in the coarse Hilbert space H b one must represent ψ b in the Hilbert space H b . Then, each states can be evaluated with their respective amplitude and the results compared. For this purpose one defines so-called embedding maps: For this to work, the boundary b must be embeddable into the boundary b , denoted as b ≺ b . Thus, the boundaries b form a partially ordered set. In case that two boundaries b and b cannot be directly related, i.e. b cannot be embedded into b , one embeds both into a common refinement b , written as b ≺ b and b ≺ b . Hence, the goal is the following: given a state ψ b in a coarse Hilbert space H b , we want to define an equivalence class of states in all finer, more complex Hilbert spaces H b in order to readily compare the associated amplitude functionals. This equivalence class of states is defined as follows: given two states For this condition to be well-defined, the embedding maps need to satisfy a consistency condition, referred to as cylindrical consistency: Essentially, it should not matter whether a state is directly embedded into a fine boundary b or via an (or any other) intermediate boundary b . Given these conditions and relations, one can (at least formally) define a continuum Hilbert space via an inductive limit: H : Beyond this formal definition, the action of embedding maps is best understood in the following way. As illustrated before, they serve the purpose of representing a coarse state in a finer Hilbert space, which can encode more complex configurations. Hence, embedding maps specify how and in which state degrees of freedom are added. Moreover, they thus define an inner product allowing us to compare states across Hilbert spaces. Since the information of the coarse state ought to be unchanged, these new degrees of freedom are added in a vacuum state prescribed by the embedding map. These concepts are familiar in the kinematical Hilbert space of loop quantum gravity [3] expressed in terms of spin network functions, where new degrees of freedom are added in the Ashtekar-Lewandowski vacuum [60,61], which describes no space. In contrast, a dual BF representation [62][63][64] constructed in the last few years adds degrees of freedom that are peaked on flat connections. However, this notion of vacuum does not imply that this is a physical vacuum. Both examples given above are kinematical vacua in 4D gravity, i.e. they do not satisfy diffeomorphism and Hamiltonian constraints.

B. Renormalization group flow of amplitudes
Once given such a choice of embedding maps, these can be readily used to compare spin foam amplitudes. Again, given two amplitudes A b and A b , a state ψ b ∈ H b and an embedding map ι b b , we compare both amplitudes: Due to the embedding map, we define an effective amplitude A b for the coarse Hilbert space H b from the fine one A b for H b . If performed for all possible states in H b , we obtain the coarse grained amplitude. Thus, embedding maps, which specify how to add degrees of freedom under refinement of states, serve as coarse graining maps for amplitudes. There, they specify how to define effective degrees of freedom. Consequently, since b can capture more information that b, embedding maps also encode how to truncate degrees of freedom.
To summarize, a class of embedding maps defines a coarse graining / renormalization group flow of amplitudes, formulated with respect to their boundary. This flow is given by the following equation: To showcase the implications for the system as a whole, it is instructive to consider the partition function of the system. For simplicity, we assume that we can write it as a collection of amplitudes A b 6 : The original partition function is given as a product of amplitudes A b assigned to building blocks with boundary Hilbert space H b . This Hilbert space is spanned by an orthonormal basis with labels j b . In the second equality, we perform a blocking of amplitudes, e.g. of 16 vertex amplitudes for hypercubic combinatorics. The degrees of freedom j b are split into two groups: one group makes up the boundary degrees of freedom j b of the blocked amplitudes, while the other are block "bulk" degrees of freedom j b and are summed over. The latter part is then summarized as the fine amplitude A b . See fig. 5. As the final step, we implement the embedding maps (or rather coarse graining maps) to derive the effective amplitude A b for the original Hilbert space.
In this context, the embedding maps serve as variable transformations and truncations, see again fig. 5. Indeed, this inclusion of embedding maps necessarily alters the partition function of the system, and we must ensure that we can still draw reliable conclusions about the original system. Thus, the combined embedding maps from neighbouring amplitudes should be close to the identity on the respective Hilbert space, in the sense that we only truncate irrelevant degrees of freedom. For example this is realized in tensor network renormalization, where these embedding maps are unitary as we explain in section V A. Essentially, with this coarse graining procedure, we achieve two goals. On the one hand, we evaluate the partition function in parts, purely from local considerations of a subset of amplitudes. This is computationally efficient and makes calculations more accessible. On the other hand, we derive an effective theory on a coarser lattice, with less degrees of freedom, from a theory defined on a finer lattice. That way, we relate two theories on two different regulators by assigning different amplitudes to different lattices. Thus, the renormalization group flow of spin foam quantum gravity is defined as a family of amplitudes assigned to a family of 2-complexes / discretisations: Before we discuss the consequences of this renormalization group flow in detail, it is important to discuss the role of the embedding / coarse graining and how they ought to be chosen.

C. Dynamical embedding maps and the physical vacuum
The embedding maps play a pivotal role in the renormalization group flow. Generically, different embedding maps result in different flows, since they determine how the effective coarse degrees of freedom depend on the fine ones. Thus, if one fixes the embedding maps a priori, this choice must be carefully checked. Instead, it is vital that these maps are directly determined from the dynamics encoded in the amplitudes themselves [30,31,40].
This reasoning is intuitive to follow, e.g. consider the Ising model. The effective degrees of freedom most suitably describing the coarse dynamics are sensitive to the temperature and significantly differ between low and high temperature. Hence, one size does not fit all: fixing embedding maps a priori can seriously distort the RG flow and give wrong results. Thus, we are convinced that dynamical embedding maps are vital for a successful coarse graining scheme. We explain how to implement this in practice in section V A on tensor network renormalization. For now, we discuss its implications and physical interpretation.
As discussed above, embedding maps are prescriptions how and in which state degrees of freedom are added under refinement. In particular they allow to relate and identify states across Hilbert spaces, defining an inner product and notion of vacuum. When these embedding maps are chosen dynamically, i.e. with respect to the dynamics encoded in the amplitude, the new degrees of freedom are added in a dynamical vacuum state.
General relativity is a totally constrained theory, i.e. in a canonical formulation its evolution is not governed by a Hamiltonian but rather by a sum of constraints, namely the diffeomorphism and Hamiltonian constraints. These constraints are generators of gauge transformations, which implies that evolution in gravity amounts to gauge transformations. This is known as the infamous problem of time [65]. For a quantum theory, the goal is to find the physical Hilbert space, i.e. the space of all states that are annihilated by all the constraints. Therefore, given an initial physical state, evolution by the constraint operators leave the physical state unchanged. Following this insight, dynamical embedding maps in a quantum gravity theory should be physical embedding maps that add degrees of freedom in the physical vacuum under refinement. Thus, such embedding maps do not add new information to the state and represent the same physical state on a finer boundary Hilbert space.
In the context of path integral formalisms of gravity, this insight is particularly intriguing: since evolution in the canonical formulation is governed by constraint operators, the path integral merely imposes the constraints, projecting out kinematical degrees of freedom and leaving physical states unchanged [30]. In short, the path integral serves as a projector onto the physical Hilbert space. Indeed, this insight is one of the original motivation behind constructing spin foam models in the attempt to define a riggin map / physical inner product for kinematical states of loop quantum gravity [3].
Following this line of thought, spin foams themselves are dynamical embedding maps. Consider a spin foam evolving an initial state to a final state, where these states are defined on different boundaries 7 . When interpreting the spin foam as a map from one Hilbert space to another, instead of as an amplitude functional, it is by definition an embedding map. Moreover, if the spin foam acts as a projector onto the physical Hilbert space, concatenated spin foams still act as a projector, implying path independence of evolution. Conversely, this is interpreted as first evolving to an intermediate state, thus cylindrical consistency conditions of embedding maps are satisfied. Additionally, the projector property implies that this evolution is independent of the choice 2-complex / discretisation, and it would mark a fixed point of the renormalization group flow: This implies that assigning the same amplitude to all discretisations gives the same results. All of the conditions mentioned above are highly non-trivial and rely on a perfect implementation of diffeomorphism symmetry in the discrete. Indeed, this assumption is hidden in the projector property of the path integral / spin foam, which implies an implementation of diffeomorphism and Hamiltonian constraints. Path / discretisation independence follow immediately and underline the strong connection of diffeomorphism symmetry and discretisation independence. In fact, this construction would be a realization of the perfect action program [38] for quantum gravity and would imply that the dynamics of quantum gravity are solved non-perturbatively and pulled back onto the discrete.
Unsurprisingly, these conditions are not met by spin foam models: spin foam amplitudes do not act like projectors [42,69] and explicitly break diffeomorphism symmetry [34,36]. Furthermore, it is unlikely that these conditions can be perfectly realized without approximations in full generality. Thus, the goal of the coarse graining scheme is to iteratively improve spin foam amplitudes in order to well approximate the ideal solution. We discuss this in the next section.

D. Lessons from the RG flow
The idea behind the coarse graining method outlined above is that it allows us to iteratively improve the amplitudes, such that the conditions mentioned above are approximately implemented. Furthermore we can tackle discretisation dependence of the theory. The expectation is that on a second order phase transition the regulator can be removed and diffeomorphism symmetry is restored. Let us explain this step by step.
Firstly, it is straightforward to recognize that the renormalization group flow addresses the question of discretisation dependence and choice of 2-complex. By deriving an effective amplitude A b from coarse graining A b , we directly relate theories on 2-complexes defined from building blocks with boundaries b and b respectively. This information is vital for any discrete theory: it states that we perform (approximately) the same calculation on b when using A b as on b when using A b . In particular, following this prescription, it does not matter whether we calculate a coarse observable on the coarse or the fine lattice. Thus, we account for the discretisation dependence of the theory and ensure at the same time that the results are reliable. Indeed, understanding this behaviour is indispensable when trying to make contact with experiments.
Understanding the lattice dependence of the theory is an important step towards determining when and how the regulator can be removed entirely. To this end, one has to study the whole coarse graining flow, that is choose an initial amplitude (e.g. given by a choice of parameters) and follow the flow until it reaches a fixed point. Within a certain approximation, e.g. restricting to finite dimensional boundary Hilbert spaces, this generically happens. These attractive fixed points frequently describe topological theories 8 , where a continuum limit can be trivially taken. However, these theories do not describe gravity in four dimensions, since they lack propagating degrees of freedom. Furthermore, these attractive fixed points denote the phases of the theory.
All initial amplitudes, e.g. all amplitudes from a certain region in parameter space, that flow under coarse graining to the same attractive fixed point lie in the same phase. These theories have the same dynamics on sufficiently coarse grained discretisations and thus lie in an universality class and share qualitative features, e.g. in expectation values of observables. An example would be the strong coupling phase in lattice gauge theory, in which one expects the Wilson loop operator to satisfy an area law. Frequently models possess multiple phases with phase transitions separating them. We are particularly interested in phase transitions of second order.
In standard lore, second order phase transitions are characterized by a diverging correlation length. This implies that degrees of freedom infinitely far away are correlated and, thus, infinitely many degrees of freedom are relevant for the dynamics. Moreover, right on the phase transition, the system is scale-invariant, i.e. physics do not change with scale. Therefore, on a second order phase transition one can take the continuum limit to arrive at a continuum theory with propagating degrees of freedom.
We expect the same to hold for second order phase transitions in spin foam models, with a slightly different interpretation: Background independent theories lack an absolute length scale. Still, the regulator allows us to define a combinatorial distance. Essentially, the idea is to define a distance between vertices of the 2-complex, by counting the number of vertices one has to pass in order to reach the other one. If they are directly connected by an edge, this distance would be one 9 . Then, on a second order phase transition, degrees of freedom that are infinitely "far" away with regard to the lattice are correlated and the combinatorial correlation "length" diverges. Furthermore, the notion of scale invariance is consequently replaced by a discretisation independence / invariance, fixed point equations for the amplitudes are satisfied and a continuum / refinement limit can be taken. The implications of constructing the theory on this fixed point must be stressed: Due to the discretisation independence, calculations can be performed in the continuum or on any discretisation, giving the same results. This exactly corresponds to the idea of perfect action, and thus solving the coarse graining flow corresponds to solving the theory on all lattices.
Nevertheless, two caveats must be observed: firstly, finding such a second order phase transition (if it exists) does not guarantee that the corresponding theory is a correct theory of quantum gravity. Secondly, if infinitely many degrees of freedom become relevant, truncated coarse graining schemes can only approximate the desired theory to a certain order, as it is the case in other renormalization schemes.
In the next section we discuss the notion of scale in more detail.

E. Background independence and the interpretation of scale
Before we continue with reviewing how to coarse grain in practice, it is crucial to discuss the notion of "scale" -or the lack thereof. As a background independent approach, one cannot assign a scale to a spin foam since one sums / superimposes all possible geometries (allowed by a certain 2-complex). Thus, we use the 2-complex itself, here formulated via the boundary of the amplitudes, to order degrees of freedom according to a relative scale. Consequently, we replace the familiar notions of ultraviolet (UV) and infrared (IR) by "fine" and "coarse" respectively. Following this perspective, we are not integrating out short scale degrees of freedom under coarse graining. Instead, we sum over finer discrete degrees of freedom and define effective coarse degrees of freedom, which encode (superpositions of) geometries of different scales.
Alternatively, one can introduce a specific scale in this coarse graining procedure via boundary states. That is, we do not consider the entire Hilbert space, but only a specific state because we are interested in studying a transition of geometries. Then, this fixed boundary states introduces a physical scale via the encoded 3D geometry, e.g. implemented in the restricted path integral formalism, see section V B.
A further comment on the coarse graining scheme is in order: here we purely formulate it in terms of the boundary discretisation and its associated boundary Hilbert space, not in terms of the bulk. From a practical perspective these questions are less important, e.g. if one assumes regular combinatorics such that the coarse graining procedure can be straightforwardly iterated; then both the boundary and bulk form totally ordered sets. Nevertheless, there is a proposal by Bahr [72] for formulating the coarse graining scheme in the bulk, essentially by defining embedding maps for 2-complexes.

V. COARSE GRAINING METHODS
For the rest of this article, let us focus on numerical methods that allow us to realize this coarse graining method in practice and review the results.

A. Tensor network renormalization methods
Tensor network methods originate in the fields of quantum information and condensed matter and aim at efficiently studying quantum many body systems. In this review we focus on tensor network renormalization methods 10 [46][47][48], a numerical algorithm for coarse graining discrete systems. To this end the partition function of the system is represented as a contraction of a tensor network. In the context of this article, a tensor T abc... is best understood as the amplitude assigned to a region. The boundary data of these amplitudes are represented as indices of the tensor, which is graphically represented as a vertex with as many legs as it has indices. The partition function is then rewritten as a contraction of tensors: Graphically, each identified and contracted index is represented by connecting the respective tensor indices. Thus, the partition function is represented by a collection of tensors connected to one another in a local fashion, a tensor network, see fig. 7.
So far, this is merely a rewriting of the original system. The goal is to locally manipulate the tensors in order to rewrite the partition function as a coarser tensor network, see again fig. 7. This may require truncations / approximations for which the error can be estimated. There exist several tensor network schemes, yet they all have a series of steps in common that we illustrate for a concrete example.
For simplicity, take a 2D quadratic tensor network. One step present in all tensor networks is an explicit summing of degrees of freedom, referred to contraction of indices. In our network we group together four tensors T and sum over their shared indices and obtain a new tensorT , which has twice as many indices, yet the network remains local, see fig. 8 We observe an immediate issue: if the original tensor had an index range of χ, called the bond dimension, the new tensor has a range of χ 2 . Thus, while we evaluate the partition function in steps, we cannot continue indefinitely without truncations / approximations. To implement those, dynamical variable transformations are derived from thẽ T via a singular value decomposition. That way, we define effective coarse degrees of freedom as functions of the fine ones. Crucially the effective degrees of freedom are derived from the dynamics encoded in the tensors. This works as follows: FIG. 7. The basic idea of tensor network renormalization: write the partition function as a contraction of a tensor network of tensors T and then locally manipulate the tensors, such that the same partition is approximated by a coarser network of effective tensors T . This defines a flow in "tensor space".
Given the new tensorT , we intend to map two indices into an effective one. To do so, we splitT in two, separating the strands a 1 , a 2 from all other variables. This is generically not possible unless the tensor factorises. To split the tensor, we first rewrite it in terms of a matrix M AB , where index A = {a 1 , a 2 } and B contains all remaining indices. On this matrix we perform a singular value decomposition: The matrices U and V are unitary and contain the left and right singular vectors of M . λ is a diagonal matrix of singular values, where λ 1 ≥ λ 2 ≥ · · · ≥ λ χ 2 ≥ 0. See right side of fig. 8. If we translate the matrix indices A and B back into the original tensor indices, we see that the singular value decomposition allows to write the tensorT as the contraction of a three-valent and a seven-valent tensor, where the summed index i labels the singular values. The three-valent tensor U encodes the desired variable transformation, translating the degrees of freedom a 1 , a 2 into an effective coarse degree of freedom / index i. This transformation is exact, since i (generically) has a range of χ 2 . Since U is a unitary matrix, we can introduce resolutions of identity U U † in the partition function (see left of fig. 9) without changing it and sum over the indices a 1 , a 2 as well as the indices c 1 , c 2 on the opposite site. Then we repeat this procedure for the remaining indices to obtain a new effective tensor T , see right side of fig. 9.
Hence, we define a new effective tensor, yet its index range is still χ 2 , and we cannot continue to iterate this procedure without truncations. The singular value decomposition allows us to implement this truncation in an optimal way. Since all singular values are positive semi-definite and ordered in size, λi λ1 indicates how significant i is with respect to the most significant one, i = 1. Indeed, we can approximate the rank χ 2 matrix M by a rank d matrix by ignoring all  The algorithm briefly sketched above is deliberately chosen to showcase that tensor network renormalization provides a concrete realization of the spin foam coarse graining scheme. Firstly, it blocks together tensors and exactly sums over their internal degrees of freedom. Secondly, the singular value decomposition provides dynamical variable transformations that fulfill the role of dynamical embedding maps, and, moreover, allow for efficient and controllable truncations. One way to check whether these approximations are justified is to gradually increase the bond dimension, i.e. the number of kept singular values, to see whether the properties of the system depend on this choice, e.g. the position of a phase transition in parameter space. In case the results converge, the approximation is sufficiently good and one can extrapolate to infinite bond dimension. However, there exist situations in which no truncation should be implemented, in particular on second order phase transitions [48]. There, one observes that more and more singular values are relevant the closer one tunes towards the phase transition, such that one would require an infinite bond dimension or equivalently infinitely fine boundary data. This is expected, since one models a highly non-local system by locally gluing amplitudes.
We would like to highlight some general advantages and disadvantages of tensor network methods, as well as modifications to the method that so far are not applied to models of quantum gravity. Compared to other numerical methods, like Monte Carlo methods, tensor network algorithms do not suffer from the sign problem; the algorithms are perfectly applicable to quantum (oscillating) amplitudes, like in spin foams. The reason is that tensor networks do not rely on (random) sampling methods for a large system, but usually focus on all possible configurations of an amplitude for one building block (tensor). However, this leads to one of their disadvantages: in order to save all configurations of a tensor, it must have a finite dimensional boundary Hilbert space. Moreover, the numerical costs, both in terms of computational time and memory usage, scale with the dimension of the boundary Hilbert space. In particular for lattice gauge theories and spin foams in higher dimensions, this requires extensive optimization to make the numerical simulations feasible.
Before reviewing results in quantum gravity (related) models, we would like to highlight a few methods from the tensor network community. One key modification is called entanglement filtering [47,48,75]. It removes entanglement between short-scale degrees of freedom, which would otherwise get promoted to larger scales and lead to unphysical fixed points in the renormalization group flow. Other modifications aim at including Monte Carlo methods into tensor network algorithms, e.g. for contracting tensor indices [76] or to sample over the probability distribution of coarse degrees given by the singular values [77].
In the following, we first review works on tensor network renormalization applied to 2D analogue models. There we focus on the introduction of symmetry preserving methods that use the symmetry of the system to label the effective degrees of freedom with the original variables. In the second part, we discuss how to apply these methods to 3D lattice gauge theories and spin foam models, which require an efficient description of the model given by so-called decorated tensor networks.

Analogue spin foam models in 2D
By 2D analogue spin foam models, sometimes also called spin net models, we mean spin systems with a global symmetry. The typical example is the Ising model (with vanishing external magnetic field) that has a global Z 2 symmetry. Typically these models are written in terms of group variables colorred g v ∈ G assigned to the vertices of the lattice, which only interact with their nearest neighbours expressed in "edge weights" ω e (g s(e) g −1 t(e) ). In order to work with finite dimensional Hilbert spaces we restrict G to be a finite group (or quantum group later on) [78]. The partition function of the system is given by: To be invariant under the global symmetry, i.e. an element h ∈ G acting on all vertices at once, these edge weights must satisfy Thus ω e are class functions and Since the function ω e are invariant under conjugation, each one can be expanded via Peter Weyl's theorem [79] into a sum over irreducible representations ρ of the character χ ρ of G: ω ρ stands for the group Fourier transform of the edge weight ω e . Performing this for all edges and expanding the characters as a trace of representation matrices, the expression factorizes over all group elements g v , such that the group integrations / summations can be performed analytically: P v denotes the Haar projector of the group G, i.e. the projector onto the invariant subspace. We suppress its many indices for clarity of the notation. After performing all group integrations / summations, the partition functions reads: Note that the indices of the Haar projectors P v are contracted with projectors on neighbouring vertices. For more details on these models and their relation to spin foam models (with finite groups) see [78]. The expansion sketched here is completely analogous to the derivation of the spin foam representation familiar from spin foam literature. Thus, while the dimensionality is lower and the dynamics simpler, the dynamical ingredients -irreducible representations ρ and projectors onto the invariant subspace P / intertwiners ι -are the same as for spin foam models. Moreover, it is expected that these 2D spin systems share statistical properties with the 4D gauge theories of the same group [80]. As a final point, these models can be related to peculiar spin foams that only possess two vertices and many edges [81].
Hence these models represent ideal test cases for applying tensor network renormalization to spin foam models and derive first hints for the RG flow of the full theory. Fortunately, the translation of the partition function into tensor network language is straightforward: the projectors P are essentially tensors, whose variables are the irreducible representations ρ e on the edges. Just the weightsω ρe need be split per edge via a squareroot. While tensor network algorithms can be readily applied, it is vital to consider the symmetries encoded in P: the irreducible representations ρ e meeting at the vertex v must satisfy the coupling rules, i.e. they must couple to the trivial representation to satisfy gauge invariance. These restrictions can be used to optimize tensor network renormalization methods in two ways: firstly, by only storing and summing over configurations allowed by the coupling rules, the memory cost and numerical cost for the singular value decomposition and index contractions can be drastically reduced. For Abelian models this is straightforward, since all representations are one-dimensional and the projector P = δ (G) ( e⊃v k e ) (modulo orientation of the edges). Thus, under splitting of a tensor, e.g. to define the variable transformations / embedding map, one defines an intermediate representation for the new effective edge satisfying the coupling rules for both tensors. This new representation will be the label of the effective degrees of freedom and thus allows us to explicitly preserve the symmetry, see fig. 10. Moreover, ordering the entries of the matrix according to the intermediate representation turns the matrix into a block diagonal form. Thus, the algorithm can be further optimized by performing a singular value decomposition individually for each block 11 . For non-Abelian models, a further comment is necessary. As sketched in equation (21), the tensor possesses "magnetic" indices m, n per edge in addition to the irreducible representation ρ. More precisely, for fixed ρ, each edge carries the vector space V ρ ⊗ V ρ * , where ρ * denotes the dual representation to ρ. Since V ρ * can be identified with the dual vector space (V ρ ) * , we label each edge with a single representation. For tensor networks these magnetic indices pose a significant challenge: if we were to include them, the size of the tensor would render the simulations unfeasible. Fortunately, the dependence on these indices is entirely encoded in representation theory of the group G and does not change under coarse graining. Essentially these indices get "pre-contracted" [81,82], which is accounted for by G recoupling symbols in the coarse graining equations. From these equations one can read off another feature of non-Abelian models: the renormalization group equations for the representations ρ and ρ * are decoupled. Thus, it is possible that under coarse graining the effective edges will carry representations (ρ, ρ ) with ρ = ρ * .
Mentioning the channels (ρ, ρ ) is a good keyword to explain the flow as well as the approximation scheme. Due to the explicit symmetry preservation, the renormalized tensors are expressed in terms of the same variables as the original theory (with a slightly more general theory space). Thus, instead of directly comparing all entries of the tensors, we study the coarse graining flow by considering the singular values per channel (ρ, ρ ). This is completely sufficient to characterize the flow and read off different fixed points. Then, in order to determine which degrees of freedom are more relevant, one needs to compare the singular values from all channels and truncate accordingly. This can then result in a higher multiplicity of the same representation labels and thus more general boundary data, which improves the accuracy of the simulation. However, in most cases [81][82][83] a simple scheme is used, where only the largest singular value per block is kept. While this is a strong simplification, it is sufficient to identify interesting phases, labeled by attractive fixed points of the flow, while keeping the simulations feasible, in particular for studies on quantum groups SU(2) k [81] and SU(2) k × SU(2) k [83].
The articles [84,85] investigate spin net models for the Abelian finite groups Z q , for so-called cut-off models. The starting point are the edge weights in the zero-temperature limit withω k = 1 ∀k ∈ {0, . . . , q}. These weights are then truncated at different levels k, which breaks the topological symmetry and it is investigated whether this symmetry is restored under coarse graining. While for low-k and high-k cut-off the high and low temperature fixed points are found respectively, there exist intermediate phases showing oscillating behaviour. In [82] tensor network methods are generalized to non-Abelian groups applied to spin nets for S 3 , the permutation group of three elements. To keep these models feasible, the coupling rules are heavily used to optimize the algorithm. The models investigated build upon a holonomy representation of spin foam models [86] and their implementation of simplicity constraints. In general they find a non-trivial phase diagram of three phases, a low temperature S 3 ordered phase, a high temperature S 3 disordered phase as well as a Z 2 ordered phase.
As a next step [81], finite non-Abelian groups are replaced by the quantum group SU(2) k with the deformation parameter q = exp( 2πi k+2 ) a root of unity 12 [87,88]. The advantage is that the integer level k defines a gauge-invariant cut-off j max = k 2 . That way it is possible to study systems with more degrees of freedom by increasing the level k, while the representation theory remains similar. Moreover, one eventually approaches full SU(2) as k → ∞. Moreover, quantum groups are physically motivated from 3D spin foam models, where they describe gravity with a non-vanishing cosmological constant [53]. The models studied in [81] are constructed from so-called intertwiner model fixed points [89], which represent topological field theories. In a nutshell, intertwiner models are "half" of a spin net model, with an edge Hilbert space of V ρ instead of V ρ ⊗ V ρ * . These models are interesting since one can directly investigate whether the two copies remain coupled or decouple under coarse graining. Indeed one finds a rich phase structure with potential second order phase transitions.
Eventually, the work [83] investigates spin net models for SU(2) k × SU(2) k that mimic the construction of 4D Riemannian spin foam models, namely the Barrett-Crane [10] and EPRL / FK model [13,14]. For the BC model several attractive fixed points are found, none of which correspond to topological BF theory. While this indicates that simplicity constraints are strongly implemented, no indications for a 2nd order phase transition are observed. This hints towards the fact that the constraints might be too strongly implemented [15]. In contrast, the EPRL / FK model shows a highly intricate flow and partially oscillating behaviour, most likely due to exciting only a few representations initially. This is a particularity of the implementation of the simplicity constraints in the Riemannian EPRL / FK model, which relate Spin(4) representations (j + , j − ) to an SU(2) representation k via j ± = 1 2 k |1 ± γ|. Note that j ± as well as k must be half integer, such that γ must be rational.
These results impressively show the potential of tensor network techniques for studying the renormalization group of spin foam quantum gravity. Moreover, they lead to the development of key optimizations and insights that are crucial for going to higher dimensional gauge systems. This is the subject of the next section.

Decorated tensor networks for lattice gauge theories and spin foams
Dimensions larger than d = 2 and lattice gauge theories pose challenges for tensor network renormalization methods. Since higher dimensional tensors carry more boundary data, the algorithm generically is more costly than its lower dimensional counterpart. For lattice gauge theories, where due to the local gauge symmetry many degrees of freedom are redundant, it is thus imperative to develop an optimal representation if one intends to cast them into a tensor network form. Moreover, these networks are generically more complex than spin systems, since several data are shared among more than two building blocks. One example are spin foams in four dimensions, where a face and the representation it carries are shared among multiple 4-simplices. A possible tensor network representation is to assign a tensor dual to each 4-simplex, yet one must introduce auxiliary tensors [90] to ensure the correct identification of shared variables. See [91] for a more extensive discussion of possible representations.
To improve on these representations, decorated tensor network algorithms are developed and introduced in [91]. The idea is to shift the perspective away from a pure tensor network representation of the system towards amplitudes with more intricate boundary data. Yet the key ideas are retained to explicitly contract bulk degrees of freedom and to dynamically define effective degrees of freedom using a singular value decomposition. Instead of tensors, represented by vertices and legs, one works with a spin foam inspired representation where amplitudes are assigned to regions, which carry boundary data, e.g. spin network data. While the assignment of amplitudes remains local, the non-local nature of gauge theories requires more complex boundary data and intricate gluing rules that cannot be cast in a simple tensor network form without introducing additional structures. The rest of the algorithm remains essentially the same: amplitudes are glued together by suitably identifying variables among them. Depending on the considered situation, some of the identified variables are not summed over and remain part of the boundary Hilbert space. On this fine amplitude one performs a singular value decomposition to derive an embedding map leading to coarse effective degrees of freedom.
The original algorithm in [91] works slightly differently by "splitting" the amplitudes explicitly. Let us briefly demonstrate the idea for the usual 2D Ising model, for which a decorated tensor network algorithm exists as well. Consider an amplitude assigned to a square given by four Ising spins σ i , i ∈ {1, . . . , 4}, A(σ 1 , . . . , σ 4 ), see fig. 11. The idea of the algorithm is alternatingly split the squares into regular triangles, such that four of these form a coarse, rotated square with a single Ising spin to sum over in their center. To do so we split the amplitude in two, separating the dependency on the spins opposite to the cut. For the singular value decomposition, we need to distinguish two sets of variables: the variables that we want to separate are encoded in the two indices of the matrix to decompose, while the shared variables will remain fixed similar to the symmetry preserving algorithm before. Hence, we perform a singular value decomposition for each configuration of shared variables, which is more efficient than a decomposing a big matrix. In a sense, this leads to a doubling of the shared variables, which is necessary for gluing them again in consecutive iterations. In our example, we get: By assigning a square root of the singular values λ i to U and V , we derive the desired amplitudes assigned to the triangles. Note that each amplitude is not just given by the configuration of three Ising spins, but also by an additional index i, assigned to the coarse edge. When combing four triangular amplitudes, the resulting amplitude for the square is given by more general boundary data, four Ising spins and four new indices A(σ 1 , . . . , σ 4 , i 1 , . . . , i 4 ), again see fig.  11. Since these indices are shared with neighbouring amplitudes, we represent them by a tensor network on the lattice dual to the squares. Thus, we have a tensor network encoding higher order corrections / more general boundary data "decorated" by the original boundary data of the system.
FIG . 11. Illustration of decorated tensor networks. Splitting an amplitude, here for four Ising spins, via a singular value decomposition generically gives rise to additional indices. These indices are understood as more general boundary data and are encoded in a tensor network dual to the lattice. This network is thus "decorated" by the remaining data.
An algorithm is developed for 3D lattice gauge theories and first applied to Abelian Z 2 lattice gauge theory in [91]. Instead of working with the original lattice, one works with the dual lattice in the strong coupling expansion. Thus the variables are irreducible representations k e of Z 2 one edge of the dual lattice and Gauss constraints (Z 2 δ-functions) on each face. The Gauss constraints can be explicitly solved to reduce the amount of data saved, and there is a freedom to choose which variables k e to gauge fix. This choice is adapted to the intended splitting.
Let us briefly sketch the coarse graining algorithm: the idea is to cut cubes in half by cutting along the diagonal of one of its faces. Four of these amplitudes are glued together to form a new (distorted) cube. This coarse graining is continued in the other directions by "rotating" the amplitude and iterating the procedure. Compared to the 2D algorithms, the splitting is slightly more complicated. In order to cut the cube along the face, it is necessary to introduce another representation along the intended "cut". This variable serves as the variable assigned to the coarse edge, and is introduced by splitting the Gauss constraint on the square face into two assigned to the two triangles. Then, the variables are gauge fixed such that equally distributed the remaining degrees of freedom are equally distributed among the split amplitudes, and some shared by both. As explained above, the shared variables will be kept fixed during the singular value decomposition and label the new amplitudes. Also, again more general boundary data arise in the form of a decorated tensor network due to this splitting, which will be assigned to edges and faces, see [91] for more details.
In its lowest order approximation, i.e. when truncating all tensor indices, the algorithm reproduces the phase diagram of Z 2 lattice gauge theory with a strong and weak coupling phase, whereas the critical coupling is found within an error of a few percent compared to Monte Carlo simulations [92]. These results are improved by keeping more degrees of freedom after the singular value decomposition, yet the computational costs grow quickly 13 . As one of the first tensor network algorithms applied to 3D lattice gauge theories it already shows promising qualitative results.
In [94] this algorithm is generalized to non-Abelian symmetry groups and applied to S 3 lattice gauge theory. While the basic idea and principle of the algorithm remains similar, it is significantly more complicated due to the non-Abelian group. The basic steps, splitting, gluing and choice of variables, are still in place, however in order to define them transformations between the holonomy and spin network representation are necessary. For the details, we refer the reader to the extensive and thorough explanation of the technical details in [94]. Let us focus instead on the results.
As for the similar work on S 3 spin nets [82], analogous simplicity constraints for S 3 are implemented in the holonomy representation [86]. The coarse graining flow is studied and three different phases found, as in [82], that correspond to a strong coupling S 3 phase, a weak coupling S 3 phase as well as a weak coupling S 3 /Z 3 Z 2 phase. The successful generalisation to non-Abelian groups as well as the derivation of the phase diagram of the theory demonstrate the potential of decorated tensor network techniques. However, this work also revealed a short-coming of using the spin network basis for labeling the boundary Hilbert space.
As we discuss in great detail in this review, at some step of the coarse graining process one sums over fine degrees of freedom. In the spin network basis, where one assigns irreducible representations to the links and intertwiners to the nodes, this implies defining an effective vertex / intertwiner by summing over representations, see fig. 12. However, usually the coupling rules at the effective vertex are violated such that the vector associated to the vertex is not an intertwiner any more and gauge invariance is broken. This is a well-known shortcoming of the spin network basis under coarse graining [95,96] and it is also expected for lattice gauge theories 14 . This can be overcome by a different representation of the boundary Hilbert space that can accommodate Gauss constraint violations (electric charges) as well as curvature excitations (magnetic fluxes). In 3D this is accomplished by the so-called fusion basis.  (2)): we can define an effective vertex by summing over the SU(2) representations associated to the inner edges. The configuration on the left is allowed by the coupling rules, thus the effective vertex on the right is allowed. However, this configuration is forbidden by SU(2) coupling rules.
Here we will only briefly sketch the main features of the fusion basis, which arises in anyon systems [97,98], (2+1)D lattice gauge theories [99] and 3D quantum gravity [100,101]. The algebraic structures are called Drinfeld Doubles, see [99,100,[102][103][104][105] for more details. Its main feature is that it diagonalizes a set of commuting operators, so-called Ribbon operators. These Ribbon operators, which contain both a Wilson loop operator as well as a t'Hooft operator, measure both the magnetic (curvature) as well as electric (torsion) excitations. Such excitations are localized on punctures carrying the magnetic and electric excitation. In lattice gauge theory, one can imagine one puncture per plaquette of the lattice. Ribbon operators surrounding a single or a collection of plaquettes then measure excitations associated with the puncture or collection of punctures.
These Ribbon operators commute among each other as long as they do not intersect. Hence, two Ribbon operators that surround a single puncture each commute with each other trivially, and they also commute with the operator surrounding both punctures. A choice of such a set of commuting Ribbon operators is encoded into the fusion basis by the choice of a fusion tree. The plaquettes are the leaves of the tree, and the connectivity of the tree determines which operators / observables are diagonalized by this choice of basis, see fig. 13. Moreover, the basis states can be transformed into one another, such that one can translate states to diagonalize the observables one intends to measure. This is a crucial concept that has a notion of coarse graining built into it. Imagine two cubes glued together: in order to derive an effective building block with effective degrees of freedom, one would like fuse the punctures of subdivided faces into one. To do so, the fusion basis must be chosen such that it diagonalizes the Ribbon operator around punctures. This ensures that the expecation values agree in both original and coarse grained case. Therefore, the fusion tree can be used to encode a choice of coarse grained observables, which is crucial for the decorated tensor network algorithm based on the fusion basis. Such an algorithm is defined for quantum deformed 3D lattice gauge theories on a cubic lattice for the quantum group SU(2) k in [49]. More details about the fusion basis for SU(2) k can be found in [100]. As in previous decorated tensor network algorithms, the basic ingredient is the amplitude associated to a cube. Its boundary Hilbert space is spanned by six-puncture states, for which a fusion tree needs to be specified. Then two cubes are glued together, which results in a cuboid with four course faces, each carrying two punctures. The goal is to compute an embedding map that fuses the two punctures on a subdivided face into one effective puncture. Generically, after gluing the fusion basis of the cuboid is not suited to do so, since it does not diagonalize the Ribbon operator surrounding both punctures. Thus, the basis must be transformed by a series of tree transformations involving SU(2) k recoupling theory, see [49] for details. Once one has arrived at a fusion basis that directly fuses the punctures as desired, one performs a singular value decomposition that splits the data of two punctures from all other variables (similar to section V A). In order to label the new puncture with the usual data, one keeps fixed the data of fusion tree directly after fusing the punctures. These data label the effective puncture, while the singular value decomposition gives the weight of these data in the effective amplitude. Also from the perspective of (coarse) observables this choice is viable, since the eigenvalue of a Ribbon operator surrounding two punctures is solely determined by these data. This allows us to use tensor network renormalization to approximately compute expectation values of Ribbon operators that we explain below.
After computing the embedding maps, they are used to coarse grain the punctures in such a way that the partition function is altered as little as possible (see again section V A). Again, only one non-vanishing singular value is kept per puncture label, such that the same boundary Hilbert space is pertained. Once all pairs of punctures are coarse grained, one obtains again an amplitude of six punctures that is associated with a cuboid. To complete one coarse graining iteration, the same procedure is performed in the other spatial directions. The cube is "rotated" by reordering the punctures and the same procedure repeated in all directions to arrive at a coarse cube with six punctures.
In [49] this algorithm is applied to 3D lattice gauge theories defined for the quantum group SU(2) k . As for quantum group analogue spin foam models [81,83], the quantum group introduces a gauge-invariant cut-off on the irreducible representations j max = k 2 . Thus, the boundary Hilbert spaces are finite dimensional and it is possible to study larger "groups" by increasing the level k (and approach SU(2) in the limit k → ∞). The lattice gauge theory is modeled via a Heat kernel action for SU(2) k parametrised by a gauge coupling parameter g. Lastly, in the initial state each puncture only carries magnetic excitations as it usually is the case in lattice gauge theory.
Let us summarize a few of the main results: at each level k there are two phases separated by a phase transition given by a critical coupling g c . For g < g c , the system flows to the weak coupling fixed point g = 0 and is thus characeterized as the deconfined phase. Conversely for g > g c it flows to strong coupling g → ∞, which describes the confined phase. The position of the critical coupling g c depends on the level k and decreases apparently linearly for small k. This tentatively suggests that for SU(2), i.e. the limit k → ∞, g c → 0 such that only the confining phase exists. Additionally, the fusion basis permits to track the appearance of electric excitations that get excited under coarse graining, even though the initial state had no electric charges. While they do not appear to be vital for the dynamics, e.g. the position of the phase transition is barely affected if electric charges are completely truncated, including electric charges is important for the behaviour of the coarse graining flow as they serve as (non-dynamical) disentangling maps. See [49] for more details.
The final result we would like to mention is the expectation value of observables, here of Ribbon / Wilson loop operators. Since the fusion basis diagonalizes Ribbon operators, it is straightforward to approximately compute the expectation value of coarse Ribbon operators, i.e. Ribbon operators surrounding a larger number of plaquettes. In lattice simulations one usually has to simulate the entire system in order to measure coarse observables. Here, we first coarse grain the amplitude to arrive at an effective amplitude for the coarse cube, for which we measure the coarse Ribbon operator around the coarse plaquette. Thus, we first coarse grain / integrate out the fine degrees of freedom and account for them (with some truncations) in the effective amplitude, for which we then calculate the expectation values of the operators. Using this method, we derive different scaling behaviours of the expectation value with the enclosed area of the plaquette, in particular we recover the area law of the Wilson loop in the confined phase.

B. Restricted, semi-classical path integrals
Despite the tremendous progress in developing tensor network methods for spin foam models and lattice gauge theories, applying them directly to spin foam models of 4D quantum gravity (either Riemannian or Lorentzian) is still out of reach, in particular for a continuous symmetry group. An attempt to make the 4D Riemannian spin foam models accessible is to study simpler models that represent a subset of the full gravitational path integral. These simplifications include restricting the degrees of freedom to specific intertwiners and representations as well as using asymptotic expansions of spin foam amplitudes valid only for large representations. Let us explain these assumptions in more detail.
Intertwiners, which determine the shape of dual 3D building blocks, can be expressed in terms of Perelomov coherent states / Livine-Speziale intertwiners [18,106]: to each face of the intertwiner one assigns an SU(2) coherent state |j, n , where j labels the irreducible representation and n is a vector on S 2 . This vector is a maximum weight state diagonalizing the angular momentum operator J n = n · J in n direction. Given these states, the coherent interwiner is given by: Each coherent state ∼ |j i , n i represents a face with area j i (j i + 1) peaked on a normal vector pointing in direction n i . The tensor product of these coherent states represents a 3D quantum building block sharply peaked on a classical geometry (if it exists) with areas and outward pointing normal encoded in the labels j i and n i . The group integration (with Haar measure dg) defines an intertwiner. Note that the spin foam partition function itself can also be expressed in terms of coherent states and an integral over the labels of coherent states, e.g. see the review [1] for more details. These coherent states play an important role in deriving semi-classical expressions for spin foam (vertex) amplitudes. When computing the vertex amplitude as a contraction of coherent intertwiners, these can be rewritten as several group integrations of contracted SU(2) coherent states. The latter part is then exponentiated and the group integration performed via a stationary phase approximation: f ⊃e j ab , − n ba |g −1 b g a |j ab , n ab =: dg e e f ⊃e 2j ab ln − n ba |g −1 b ga| n ab .
Since the stationary phase approximation is only valid when the argument in the exponential is highly oscillating, all representation j f must be large. Hence this expansion is often called the large-j limit. For single vertex amplitudes it is shown that the "action" in the exponential evaluated on stationary and critical points is given by the Regge action of the building block dual to the vertex [22,23]. Given these familiar results from spin foam literature, the idea is to restrict the spin foam partition function of the EPRL / FK model to specific coherent intertwiners (and representations) and use only the amplitudes derived in the large-j limit. That way the system depends on significantly fewer variables and the spin foam amplitudes, in particular the vertex amplitude, can be expressed in terms of closed formulas of the representations. Additionally, in the large-j limit the sum over representations can be approximated by an integral. The motivation is to employ numerical integration techniques, e.g. the Cuba package [107]. From now on, we only consider and discuss models defined on a 2-complex with hypercubic combinatorics which makes iterating the coarse graining steps straightforward.
So far, two models of restricted spin foams are defined that are also studied under coarse graining. The first one are so-called quantum cuboids [36], where the intertwiners are sharply peaked on a classical cuboid geometry. Opposite faces of the intertwiner carry the same representation and their normals are anti-parallel. Moreover, the outward pointing normal of a face is orthogonal to all normals of the adjacent faces, see the left part of fig. 15. Indeed these are severe restrictions, in particular the requirement that opposite faces in each intertwiner carry the same representation translate through the entire lattice. The asymptotic expansion of the vertex amplitudes depends again on the Regge action, which generically vanishes for cuboid configurations. For larger complexes this implies that the flat cuboid building blocks are glued in a flat way. Thus, this model describes a superposition of flat discrete space-times of different distribution of sizes across its building blocks. While this is by no means a realistic model of quantum gravity, it captures an Abelian subgroup of diffeomorphisms corresponding to shifts of entire hypersurfaces in the lattice 15 . Notably, this spin foam model is not invariant under these transformations [36]. Due to its simplicity it does not have any free parameters, thus an additional parameter α is introduced in the face amplitude of the model, (2j f ) α , which can be understood as a modification of the path integral measure. This exponent simply emphasizes small or large representations / face areas in the partition function, and is motivated by a discussion in the community on the right choice of face amplitude [108]. A physically more interesting model is based on so-called frusta [109]. A frustum is a higher-dimensional analogue of a trapezoid. In 4D, it consists of two cubes at its top and bottom, potentially of different size, which are connected by six 3D frusta, see the right of fig. 15. That way, a hyperfrustum can describe the evolution from one spatial cube to a larger / smaller spatial cube. Thus the idea is to restrict the intertwiners to be cube / frusta shaped in order to study the expanding / contracting cubulated 3D spatial slices. The frustum intertwiner is then given by three representations, j i and j f correspond to the initial and final area of the initial and final cubes respectively and k gives the side-face area, which also determines the "opening angle" φ of the frustum 16 . Crucially, in contrast to the cuboid model, the Regge action associated to a hyperfrustum in the asymptotic expansion no longer vanishes and one obtains the familiar cosine formula [109]: where S R denotes the Regge action of the hyperfrustum, V its volume, G Netwons' constant, Λ the cosmological constant and γ the Barbero-Immirzi parameter. Thus, this model captures several important generalizations compared to the cuboid model. As signified by the non-vanishing Regge action, frusta configurations allow for curvature to appear. Moreover, more parameters play a role in the dynamics: Newton's constant G enters here as providing an explicit scale to the representations / areas on the boundary that serve as initial and final states. The cosmological constant Λ is added in [110] (analogous to [111]) by deforming the vertex amplitude. The parameter α remains as in the cuboid model 17 .
After this basic introduction of these models, let us discuss the simplified coarse graining scheme and results.

Coarse graining setup and results
While in spirit the coarse graining setup is similar to the general method outlined in section IV, there are several noteworthy difference and assumptions being made. Firstly, the embedding maps are chosen on geometric grounds instead of determining from the dynamics / the amplitudes. The intuitive idea is, e.g. in case of the hypercuboids, that a coarse hypercuboid arises as a superposition of fine hypercuboids consistent with the coarse hypercuboid geometry. Secondly, the coarse graining flow is computed for one fixed coarse boundary state, not the whole coarse boundary Hilbert space. Thus, in the path integral context, this coarse graining flow is performed for a fixed transition. The third and final assumption includes a projection back onto the original form of the amplitude, such that the flow is formulated as a flow in parameter space of the theory. This projection is defined by comparing expectation values of observables in the coarse and fine calculation. In a sense, the logic of the consistent boundary formulation is inverted: the RG flow assigns a family of amplitudes to different lattices such that expectation values of observables agree on all lattices. Instead, we derive the RG flow by identifying theories / parameters across lattices for which the expectation values of some (sensitive) observables agree. In case of the cuboids, which are just given by the parameter α, this would read: which defines the flow α → α from fine b to coarse b. Note that all these three assumptions are strong simplifications that need to be lifted to verify their validity. Let us first discuss the setup for coarse graining in the quantum cuboid model reported in [112]. Consider two hypercuboids glued together along a common 3D cuboid. The total geometry of both hypercuboids is fixed in the coarse boundary state, i.e. the total area of each coarse face is fixed, yet the distribution of 4-volume among the two cuboids fluctuates. Obviously, the expectation value of the volume of a single hypercuboid is always exactly half of the total volume, yet its variance depends sensitively on the parameter α. Thus, for fixed coarse boundary state, one studies the variance of a single coarse in the coarse and in the fine case, where each of the hypercuboids is subdivided into 16. The geometric embedding map is prescribed such that the fine areas sum up to the coarse area. This setup and observable is particularly interesting since it is closely connected to the Abelian subgroup of diffeomorphisms that can be represented in the quantum cuboid model.
With this setup, the variance of the 4-volume is computed in both the coarse and fine case for various α / α respectively. In both cases the observable is monotonously decreasing and both curves intersect once in the value α * . This particular value of α defines a fixed point of the renormalization group flow α → α, which is repulsive, i.e. α > α for α > α * and α < α for α < α . Thus, the fixed point also separates two phases, which are dominated by different configurations. For α < α * , small representations j are preferred, such that subdivided faces contain one large area and several small ones. In contrast for α > α * , the configuration dominates in which a face is equally subdivided since then all spins are as large as possible. Remarkably, the value α * is close to the one at which diffeomorphism symmetry is almost restored [36]. This result together with the repulsive behaviour at the phase transition indicate that this transition might be of second order, and that on it the subgroup of diffeomorphisms might be restored.
The same calculation is repeated for different coarse boundary states in [113] and results in the same qualitative behaviour, yet the position of the fixed point changes. Thus, we do not include the exact value. This result sheds a light on the possible interpretation of these results. Since the coarse boundary state is kept fixed, this coarse graining derives a family of amplitudes on a family of lattices for this specific transition. Therefore it contains the information whether and for which parameters the regulator / the lattice can be removed and the results are consistent (within the given approximations and truncations). Note however that this is a weaker condition than the coarse graining flow defined in section IV, which refers to all transitions / boundary states. In a sense, the fixed state becomes part of the observable for which the coarse graining flow is defined. That way it provides first insights of coarse graining flow in a truncated theory space.
A similar analysis of the coarse graining flow is performed for frusta spin foams in [114] with a slighly different setup. Here the boundary is made up of two parts, an initial and a final 3D spatial cubulations each prescribed by a single representation j i and j f , which are chosen to be equal. Again, the goal is to compute expectation values of observables in a fine and a coarse setting and define a renormalization group flow in parameter space (α, G, Λ) such that these observables agree. In [114] a few different setups are examined, here we just discuss the main result of the RG flow in three dimensional parameter space.
Due to the high symmetry of frusta configurations, the lattice is prescribed by spatial and temporal subdivisions. The former fixes the fineness of the spatial cubulation while the latter determines the number of time steps. The coarse lattice then has two spatial subdivisions, i.e. 4 3 spatial cubes, and one intermediate time step. The fine lattice has one more spatial and temporal subdivision, i.e. 8 3 cubes and two intermediate time steps. The boundary states in both settings are straightforwardly related by requiring that the total 3D spatial volume encoded in initial / final state agrees in both settings. Moreover, the total "height" / "time" is fixed in both settings as well as each time step is chosen to be equal. That way, only the intermediate spatial volume is integrated over, while the side panels are fixed, which greatly reduces the numerical cost.
In order to derive a renormalization group flow in a parameter space with three parameters, one must consider at least three observables. In [114], the 3D spatial volume of the intermediate slice, its variance as well as the total 4D volume are considered. The expectation values for all these observables are computed in a range of all parameters G, Λ and α and compared for both settings. Then a coarse graining flow is derived by matching theories with the smallest relative error of observables 18 . Under this premise, indications for a fixed point around α * ≈ 0.677, G * ≈ 0.037 and Λ * ≈ 0.08 are found. While the exact numerical values are less relevant and most likely subject to change for different boundary states, qualitatively the numerics indicate that this fixed point has one repulsive and two attractive directions. As for the cuboid case, the repulsive direction appears to be (mostly) related to the parameter α, while G and Λ seem to be the attractive directions. In standard lore, this would imply that both G and Λ are irrelevant couplings and fixed by the RG trajectory.
At first sight, this result appears to be at odds with results in asymptotic safety [115], where both G and Λ correspond to free parameters / relevant directions. However note that this setting here is significantly different: due to considering only a specific transition, a scale is introduced into the system, which is not changed by the coarse graining flow. Thus, in contrast to asymptotic safety where one derives theories at different scales, this coarse graining flow teaches us whether and for which parameters the regulator / lattice might be removed. In this sense, the fixed point gives the correct discretisation independent amplitude (given the introduced approximations) for this specific transition. So G * and Λ * mark the correct parameters for one specific transition and are thus irrelevant in this flow, yet they might correspond to relevant directions when different scales are related.

Numerical methods
A short comment on the numerical methods is in order. In the semi-classical, restricted spin foam models, Monte-Carlo and numerical integration techniques are used [107]. This works particularly well for the quantum cuboid model, where the action vanishes and the amplitudes do not show an oscillating behaviour. Nevertheless, also for the frusta models, which feature oscillating amplitudes, can be explored with these methods with slower convergence. In general, convergence slows down for higher dimensional integrals and larger discretisations, such that this method appears to be feasible for systems with a few building blocks and symmetries that reduce the amount of degrees of freedom.
C. Semi-classical continuum limit Before concluding this review, we would like to briefly discuss the semi-classical continuum limit approach [116,117], since it aims at defining a flow across a family of triangulations and might at first sight be similar to the restricted path integral method.
This approach discussed in the papers [116,117] aims at defining a semi-classical continuum limit for spin foam models in the following sense. As mentioned before, it is a well-established result that one obtains (area) Regge calculus in the asymptotic expansion of spin foam vertex amplitues [21][22][23][24][25]. Often this is called a semi-classical limit, by scaling all representations j → λj by a parameter λ. In [117] a Gaussian weight is introduced into this semi-classical limit that suppresses non-length-Regge like geometries, i.e. geometries prescribed by areas that do not correspond to triangulation given by edge lengths. The parameter for this Gaussian weight is δ, where δ → 0 removes this weight. These systems are studied in the regime λ δ −1 1, such that the semi-classical formula is valid and higher curvature corrections (in the deficit angle) are suppressed.
Essentially the idea is to define a continuum limit as in Regge calculus: for a sequence of triangulations K N of the same manifold continuum general relativity is restored if all lengths and deficit angles converge to zero as N → ∞. In [117] it is argued that this is achieved for particular scaling relations for λ, δ and µ, where µ rescales the Planck length. Thus, one defines a flow across triangulations in this parameter space, where in the limit N → ∞, both areas and deficit angles converge to zero. Invoking the continuum limit of Regge calculus, it is argued that general relativity is obtained in this limit.
The existence of such a regime, where one obtains general relativity as the continuum limit of Regge calculus, would be intriguing and it is suggestive to think it should exist, given the close relation of spin foam models and Regge calculus. However, there are several points that must be considered before this can be confirmed: firstly, the assumptions and modifications made that must be carefully cross-checked. Most strikingly, a term that suppresses non-Regge like geometries is not present in spin foam models and one might argue that such a role ought to be already implemented in the simplicty constraints. Secondly, the conditions under which the formulas are valid are highly specific and it must be validated whether these are satisfied in generic situations. Finally, the defined flow of parameters is not dynamical, in the sense that it is not derived by relating dynamics across different triangulations. Thus, it is not clear whether this continuum limit gives well-defined continuum dynamics.

VI. OUTLOOK: TOWARDS RENORMALIZATION IN 4D
In this article we provide a detailed review of coarse graining in spin foams at the conceptual and practical level. Attentive readers notice that these methods have not been applied yet to the full 4D theory, e.g. the EPRL / FK in the Riemannian or Lorentzian setting, and it is currently out of reach. In this outlook, we would like to discuss the open issues and questions that need to be addressed.
A first point, which is relevant for all calculations performed in spin foam quantum gravity, is the computability of spin foam amplitudes, more precisely the vertex amplitude. As the amplitude associated to a 4D building block, it is the centerpiece of the theory and the most intricate to compute. Analytical formulas are known for the asymptotic expansion of the amplitude [21][22][23][24][25], where the boundary data is given by coherent states peaked on classical discrete geometries. However, these results are not valid for small representations, the quantum regime of the theory. To compute the amplitude in this regime requires numerical techniques, e.g. by explicitly contracting intertwiners to obtain the vertex amplitude. Significant progress was made in recent years for the Lorentzian EPRL model in [28,29,45] using the results form [118]. Nevertheless, these calculations require significant numerical resources, which makes it difficult to explore systems with multiple vertex amplitudes. Two ideas might be helpful in exploring larger 2-complexes: Firstly, storing computed vertex amplitudes, e.g. for an orthonromal basis of intertwiners, in an open-data database such as the "Encyclopedia of Quantum Geometries" 19 would make them accessible to interested researchers and avoid computing the same amplitudes multiple times. The second idea relies on the fact that the asymptotic formula well approximates the vertex amplitude for fairly small representations, j ∼ 10 for a 4-simplex in the Riemannian EPRL model [119]. Exploiting this fact could lead to an efficient hybrid algorithm, similar to the idea used in loop quantum cosmology [120], that only uses the costly to compute quantum amplitude in case the asymptotic formula is not accurate.
An alternative route towards studying spin foams with multiple simplices lies in defining simplified models. In section V B we review one example for such models, namely restricted spin foam models. Instead of exploring the full spin foam path integral, only a subset of configurations is explored using the asymptotic formula. Thus the number of degress of freedom is drastically reduced and the issue of exactly computing the vertex amplitude is circumvented, which makes it possible to renormalize these models. Clearly, as next steps these restrictions need to be lifted in order to explore more of the dynamics of the theory. This could either be by allowing more configurations in the path integral, see e.g. [121], or by going beyond the asymptotic formula and including the full vertex amplitude. Recently, another simplified model has been constructed in a similar direction [122]. Again, the asymptotic formula of the vertex amplitude is invoked to define a simplified vertex amplitude. Special emphasis is given to an implementation of simplicity constraints akin to spin foam models as weak conditions on 3D dihedral angles, which might give new insights into spin foam dynamics for large 2-complexes. Note that this model does not restrict the allowed configurations in contrast to the restricted models discussed in this review.
The most holistic approach to coarse graining spin foam models, tensor network renormalization discussed in detail in section V A, faces two main challenges when going to 4D. One is the increased complexity of the amplitude which results both in larger memory cost as well as computational time. Related to this is the second issue, how to define a tensor network algorithm for systems with infinitely many configurations or continuous variables like the 4D spin foam models defined for Lie groups. A solution to the former challenge might lie in defining a representation of the model suited for renormalization, similar to the fusion basis in 3D [49,100]. Alternative formulations of 4D models are investigated e.g. in [123,124]. Using observables might again serve as a guiding principle to find such representations. The second issue might be tackled in a similar direction, where a rewriting of the model might lead to an efficient tensor network description. One such example is [125], where the renormalization group flow of φ 4 scalar field theory is accessible for tensor network methods by performing a simple transformation.
Another important research direction, on which renormalization and coarse graining can shed a new light, is matter coupling in spin foam quantum gravity. Since spin foams are a purely gravitational theory, matter degrees of freedom must be added in order to adequately describe the universe. Different ways to couple matter to spin foams exist in the literature [126][127][128][129][130][131], yet the intriguing dynamics of the coupled matter gravity system are hardly explored. Applying a coarse graining scheme to the combined system allows us to renormalize matter and gravitational degrees of freedom at the same time, uncovering the phase diagram of the system. This idea is realized for a simplified toy model in [93]. Without a question, a system consisting of both spin foam and matter degrees of freedom is more difficult to study than the former alone. Nevertheless, adding matter to simplified models might be accessible and lead towards intriguing new features and insights, e.g. it would be interesting to see how the matter sector influences the quantum gravitational theory as in [132].
Beyond the methods discussed in the review, ideas from other fields and approaches to quantum gravity might help us advance coarse graining in spin foam models to 4D. These might be novel numerical techniques, like deep learning, or well-established ones like Monte Carlo methods, which might be efficiently applicable in certain settings.

CONFLICT OF INTEREST STATEMENT
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.