Renormalization of Group Field Theories for Quantum Gravity: New Computations and Some Suggestions

We discuss motivation and goals of renormalization analyses of group field theory models of simplicial 4d quantum gravity, and review briefly the status of this research area. We present some new computations of perturbative Group field theories amplitudes, concerning in particular their scaling behavior, and the numerical techniques employed to obtain them. Finally, we suggest a number of research directions for further progress.

Beside a quadratic local term (defined by an integral kernel convoluted with the two fields), the GFT action is determined by interactions (just one in the above example) that possess a characteristic 'combinatorial non-locality'in that the interaction kernels pair non-locally the field arguments (d variables being contributed by each field entering the given interaction term). Interaction kernels of order n can be associated with possible ways of gluing together n (d − 1)-simplices to form (the boundary of) a d-dimensional cell. The specific combinatorial patterns (i.e., the specific cells being associated to each fundamental interaction) and precise form of the interaction and kinetic kernels are part of the definition of each particular GFT model. However, from this generic aspect of the formalism follows one key fact: GFT Feynman diagrams Γ generated by the perturbative expansion of the GFT partition function, A Γ obtained gluing interaction vertices (d-cells) along their (d − 1)-faces, are dual to d-dimensional cellular complexes, of arbitrary topology (since a priori there is no restriction on the allowed gluing). The Feynman amplitudes assign a probability amplitude for each such cellular complex, seen as an elementary interaction process of the fundamental GFT quanta. The above can be taken as a sketchy definition of the TGFT formalism, but it is of course the specification of particular models which gives tentative physical meaning to it and this meaning will therefore change in different contexts. In particular, when the physical interpretation of a given model is grounded in its perturbative expansion, it will affect what we expect about the properties of the Feynman amplitudes A Γ . In this contribution, we focus on GFT models for quantum gravity and in particular on the class of models closely related to simplicial gravity path integrals, spin foam models and loop quantum gravity. For this class of models we discuss the key features of the Feynman amplitudes in the next subsection. In particular, beyond their interpretation, we will discuss at length whether we should expect them to be finite or divergent when dwelling into the issue of renormalization of the corresponding quantum gravity models.
Before doing so, let us spend some words to clarify the choice of a statistical form for the GFT partition function. The foundations of the GFT formalism have received some attention only recently, and much remains to be understood. Both statistical and complex weights for GFT fields summed over in the definition of the partition function can be found in the literature, but in fact most of the literature until recently has been focusing only on the perturbative expansion of the models, where the non-perturbative definition is less directly relevant, thus avoiding the need to investigate it. In particular, for the class of models to be analyzed in the next subsection, the interpretation in terms of quantum sums over discrete geometries (and topologies) is guaranteed by the form of the Feynman amplitudes (which are either complex or real with oscillating behavior, in line with their interpretation as defining a quantum discrete path integral for gravity on a lattice, as we are going to discuss). As for the non-perturbative definition of the GFT dynamics, instead, we have less guidelines, especially for the general definition of the formalism, before considering specific models. In the context of tensor models and tensorial group field theories treated as a generalization of usual local QFT framework, all non-perturbative analyses have simply assumed from the start a statistical definition, taken to be primary and in no special need of further justification; in particular, no attempt to understand it as the counterpart of an operatorial definition of the same model has been made. In non-quantum gravity applications, this is simply a choice of a 'classical'statistical field theory context, usually adopted for purely mathematical reasons (better chance of making sense of the path integral).
For quantum gravity applications, one may be interested also in a better conceptual foundation of the formalism. In this context, a complex weight involving the GFT action as a pure phase is of course a possibility, and maybe more in line with the intuitive idea of a quantum gravity path integral (even though the quantum gravity interpretation of the GFT field itself is not straightforward, while as we remarked already, this interpretation be consistently associated to the discrete structures appearing in the GFT perturbative expansion). However, we should remind ourselves that these quantum gravity models should not be expected to encode any global unitary evolution, as it is in fact true for any fundamental quantum gravity dynamics, and this removes one strong motivation to insist on this formulation of the path integral. Nor we have a complete derivation of the GFT path integral from a canonical quantum gravity dynamics or, alternatively, from some formal field theoretic gravitational path integral, that could dictate one choice over another. To date, the only tentative derivation of a GFT partition function 'from first principles', was given in ref. 4 where it was seen to arise from a quantum statistical definition of equilibrium for a system of quantized simplices (indeed, the basic quanta of such GFT models), under a requirement of maximization of entropy and a choice of macroscopic conditions to be imposed on average (this choice concurs to the specification of the resulting GFT action). If among these constraints one includes some appropriate counterpart of the Hamiltonian constraint of canonical gravity (adapted to the discrete setting), then the GFT partition function can be seen as a sort of grandcanonical partition function relaxing the imposition of such constraint (in the sense that configurations satisfying the constraint are assigned greater weight, but fluctuations off the constraint surface are allowed) as compared with a 'microcanonical'ensemble in which only solutions of the constraint are allowed. From a canonical quantum gravity perspective the latter, more restricted case would correspond to a definition of the physical inner product between quantum gravity states (with appropriate insertions of observables inside the GFT partition function), while in general the GFT formalism deals then with a broader class of quantum amplitudes. This scenario was also anticipated (more formally) in ref. 5, discussing the relation between GFT and canonical loop quantum gravity. This perspective also resonates (with many details still to be clarified, though) with the presentation of GFT from a quantum gravity perspective in ref. 6, where the 'tree level'(thus dominant, from the perturbative point of view) GFT amplitudes were suggested to define the physical inner product of canonical quantum gravity, while the remaining GFT configurations were associated with topology changing processes (off-shell, from the canonical quantum gravity perspective).

Group Field Theories, Spin Foams and Other Quantum Gravity Formalisms
We recognize in this brief outline the straightforward generalization of how 2d surfaces are generated in the perturbative expansion of random matrix models. Indeed, GFTs can be seen as group-theoretic enrichment of random tensor models [7][8][9], to which they reduce if the Lie group domain is replaced by any finite set of N elements. The Feynman amplitudes become purely combinatorial, but the type of diagrams remains the same. Seen as tensors, GFT fields admit a natural action of unitary (and orthogonal) groups on their arguments. If one requires GFT interactions to be invariant under such unitary transformations, they can be fully classified, and we speak of tensorial GFT models. Most of the literature on GFT renormalization [10,11] concerns these tensorial GFT models. More generally, focusing on tensorial aspects of GFTs allows to gain a greater control over the combinatorial structures of their states, diagrams and amplitudes and many of the results obtained in the simpler context of tensor models apply also to GFTs: the use of colors to encode the topology of Feynman diagrams, the large-N expansion, double scaling limits, universality results etc. The first two, in particular, are crucial for GFT renormalization.
In this contribution, we focus on GFT models which are 'quantum geometric': their fundamental quanta are quantized tetrahedra with a quantum geometry encoded in group-theoretic data. More precisely, the classical phase space of a single Lorentzian tetrahedron in 4d is chosen to be the cotangent bundle of 4 copies SL(2, C), reduced by additional 'geometrician'constraints, and in turn this can be mapped, under the same constraints, to the cotangent bundle of 4 copies of SU (2). At this classical level, this map amounts simply to a change of variables between two alternative parametrizations of the same classical geometry of an individual tetrahedron. We will give more details on the simplicial geometry in the next subsection (see also ref. 12). Appropriate gluing of five geometric tetrahedra on the boundary of a combinatorial 4-simplex can then be shown to provide a geometric characterization of the 4-simplex too, and of the whole simplicial complex obtained gluing geometric 4-simplices together. The same construction in the Riemannian case uses Spin(4) instead of SL(2, C). From the choice of classical phase space follows a choice of Hilbert space for individual quanta of the GFT model, given in one representation by L 2 (G 4 ) reduced by the quantum counterpart of the geometrician constraints, where G is one of the chosen groups mentioned above. More precisely, the natural Hilbert space for a single tetrahedron in this class of models would be L 2 ((SL(2, C) 4 ) restricted by the geometrician conditions dictated by the underlying classical simplicial geometric understanding of the GFT quanta (and amplitudes); this Hilbert space can be mapped, however, to the Hilbert space L(SU(2) 4 ), with the geometrician constraints encoded, in a model-dependent manner, in the definition of the map. In this latter case, the covariance properties of states and amplitudes under the action of the Lorentz group as well as some of their geometric features become 'hidden'in the form of the kernels defining the GFT action or in the corresponding spin foam amplitudes, while boundary data (and GFT fields) only depend on SU(2) data. Spin foam and GFT models defined using one or the other choice of fundamental Hilbert space are, in general, not equivalent, but the precise relation depends on the properties of the map being used (for example, its being isometric or not), but they would be equally justified from a simplicial geometric point of view. A definition of such map for the EPRL model has been given in ref. 13, together with an analysis of its properties, and a generalized definition of such map and analysis of its properties, valid for the whole class of models we deal with here, can be found in ref. 14. The complete GFT Hilbert space is the corresponding Fock space built on this single-quantum Hilbert space. We will give a few more details in the next subsection.
For this class of geometric models, the GFT formalism benefits from direct links to other modern quantum gravity approaches, which can, viceversa, benefit from GFT tools and results.
First, of all, when SU(2) is used, the Hilbert space of a single GFT quantum is the same as that of a loop quantum gravity (LQG) [15] 4-valent spin network vertex. Generic GFT states, organized in a Fock space, will be populated by many such vertices and they will include, in particular, states corresponding to spin networks associated to closed graphs and gauge-invariant cylindrical functions for the same graphs. In fact, the LQG Hilbert space associated to any graph can be shown to be faithfully embedded in the GFT Fock space. The theories however differ in the way these graph-based Hilbert spaces are related, more precisely, in the scalar products between states associated to different graphs. Still, the correspondence, which can be extended to observables and quantum dynamics, allows to see GFTs as a 2nd quantized counterpart of LQG [5].
Next. for this class of models, the GFT Feynman amplitudes take the form of (non-commutative) simplicial gravity path integrals [16,17], when written in (non-commutative) Lie algebra variables, which encode the discrete metric. The group variables, on the other hand, are understood as encoding the discrete gravity connection. They correspond indeed to discretizations of a classical formulation of gravity as a topological BF theory with added geometrician constraints, on the simplicial complex dual to the GFT Feynman diagrams. The specific way in which the BF action is discretized depend on the quantization map applied to Lie algebra variables, and different models correspond to different strategies for the imposition of the constraints and path integral measures.
In fact, when the same Feynman amplitudes are recast as functions of group representations, using Peter-Weyl or Plancherel decomposition, they take the form of spin foam models [18]. Spin foam models have been introduced as a covariant language for computing spin network dynamics, so they can be understood as a covariant counterpart of canonical loop quantum gravity. A second perspective is to see spin foam amplitudes as a purely algebraic version of lattice gravity path integrals, or state sum models. In GFT, they arise as Feynman amplitudes. The correspondence is generic: for any given set of spin foam amplitudes associated to simplicial complexes (and admitting a local decomposition with respect to the complex), one can find a GFT action such that the perturbative expansion of the quantum partition function will produce the given amplitudes as Feynman amplitudes (and viceversa, any GFT action corresponds to a set of spin foam amplitudes). A complete definition of a spin foam model requires a prescription for the amplitudes to be associated to all possible cellular complexes (in some specified class) and an organization principle for them, i.e. one way of comparing, composing or selecting them, to obtain a single number for any observable one wants to compute. The GFT embedding provide one such clear organizing principle, by summing them in a QFT perturbative expansion. In addition, it provides a whole set of QFT tools that can be applied to study their mathematical foundations as well as for extracting physics. GFT renormalization can be seen, indeed, from this spin foam perspective.

Simplicial Group Field Theory Models for 4D Quantum Gravity
The starting point for the construction of simplicial GFT (and spin foam) models of 4d quantum gravity is the quantum geometry of a single tetrahedron in 4d [19].
The quantum geometry of this basic building block, and the extended structured built from it, can be described in various parametrizations [20,21], and a number of generalizations can also be defined [22,23] and imported in the GFT framework. Classically, one can use two equivalent characterization of a tetrahedral geometry, leading immediately to an algebraic translation. First, one can start with assigning four vectors b I i ∈ R 3,1 to the four faces of the tetrahedron, forced to lie all in the same spacelike hypersurface with timelike normal V (thus satisfying b i · V 0), and thinking of them as normal to the same faces, with their modulus identified with their area, b I i A i n I i (with |n| 1). The vectors are also forced to close to form the closed boundary of the tetrahedron, i.e. i b i 0. The vectors b i , due to the constraints they satisfy, are actually elements of the vector space R 3 which can be identified with the Lie algebra su(2), after it has been endowed with the corresponding Lie bracket. The resulting space su(2) ×4 is then the space of geometries for a single tetrahedron. It can also be seen as the cotangent space of the phase space (T * SU(2)) ×4 which is then the phase space of a classical tetrahedron, purely expressed in terms of group-theoretic, algebraic data. The conjugate variables in SU(2) ×4 have the interpretation of parallel transports of a discrete connection along elementary paths from (the (bari)center of) the tetrahedron to the ((bari)center of its) boundary faces. The dual graph made of these paths becomes the graph associated to a single spin network vertex (with four outgoing 'open links'). In group representation, the corresponding Hilbert space is thus L 2 (SU(2) 4 ) (with Haar measure).
An equivalent encoding of the classical geometry of a single tetrahedron uses directly the variables of discretized topological BF theory. All geometric quantities of a single tetrahedron can be computed starting from four bivectors B IJ i ∈ ∧ 2 R 3,1 xsl(2, C) which close i B IJ i 0 and satisfy the simplicity constraints V I (*B i ) IJ 0 (* is the hedge dual), with respect to the same timelike normal vector V. The phase space of a single tetrahedron can be taken to be the cotangent bundle T * SL(2, C) 4 and the Hilbert space to be L 2 (SL(2, C) 4 ). See also ref. 12 for more details. This second construction can be indeed seen as the discrete (and then, quantum) counterpart of the formulation of continuum General Relativity as a constrained BF theory (a topological field theory) in 4 spacetime dimensions. This amounts to adding suitable constraints, called 'simplicity constraints', to the BF action, resulting in the B field of the topological theory being equivalent to a tetrad field, in such a way that the insertion of the general solution of these constraints in the BF action gives the Palatini formulation of classical continuum gravity in terms of tetrad and connection fields, in turn equivalent, at the classical level, to the metric formulation (modulo subtleties concerning degenerate geometries). The 'geometrician'constraints we discussed above correspond to the combination of the (discrete counterpart of the) simplicity constraints and the gauge invariance constraints. For more details on the continuum formulation, see the cited references. The two geometric descriptions can be mapped into each other, as we have mentioned already. The simplicity constraints can be seen also, in fact, as determining such map [13,14]. Since the Hilbert spaces indicated above admit a basis labeled by group representations, this correspondence can be seen also at that level, i.e. as specifying how the relevant representations of SL(2, C) should be decomposed in SU (2) representations, if they have to be understood as encoding the quantum geometry of a tetrahedron. Such representation labels are the variables in which spin foam amplitudes are expressed. Different spin foam (and GFT) models for 4d quantum gravity are specified (among other things) by the way the impose the simplicity constraints at the quantum level, and thus by the specific map between SL(2, C and SU(2) entering their amplitudes, if used. In the Riemannian case, which will be our focus in this contribution, all the above applies, with Spin(4) replacing SL(2, C).
A spin foam amplitude, that is a GFT Feynman amplitude written in representation variables, will be assigned to any given simplicial complex, dual to a GFT Feynman diagram. The basic building block is an assignment of a quantum amplitude to each 4-simplex, i.e., a 'vertex'of the spin foam complex given by the GFT Feynman diagram, with this amplitude function of the algebraic data associated to the five tetrahedra on its boundary. These boundary data can be written as SU(2) or SL(2, C) data, using the mentioned map, and the vertex amplitude can be written as a function of both, featuring then the coefficients of the map, if used, and the geometrician constraints, in its expression. Thus the vertex amplitude will be a function of SU(2) and SL(2, C) representations associated to the triangles of the 4-simplex (faces of the dual complex), and intertwiners of both groups associated to the tetrahedra, following Frontiers in Physics | www.frontiersin.org February 2021 | Volume 8 | Article 552354 the imposition of the closure conditions (equivalent to gauge invariance with respect to both groups). The data not used as boundary data are then summed over independently in each vertex amplitude. The spin foam amplitude associated to the whole simplicial complex can then be obtained by gluing together the amplitudes associated to its 4-simplices, with the gluing amounting to matching first and then tracing over the data associated to the tetrahedron shared by each pair of 4simplices, possibly weighted by an additional gluing kernel. In the GFT context, the vertex amplitude and the gluing kernel are nothing else than the interaction kernel and the propagator (inverse of the kinetic kernel) defining the GFT action. The correspondence between GFT amplitudes and spin foam models, which could be motivated and defined independently, is thus very general [24].
One last comment about the discrete geometry of these models. The construction sketched above, at the classical level, leads to a full characterization of the discrete geometry of the 4d simplicial complex (to which the spin foam amplitudes are associated), equivalent to the more standard characterization in terms of edge lengths, as used in Regge calculus, even though it uses a different set of classical variables (it corresponds, indeed, to a formulation of classical simplicial geometry in terms of the discrete counterpart of the variables of BF theory, suitably constrained, or to so-called 'areaangle'Regge calculus [25]). The translation of the same characterization at the quantum level, and in particular the correct imposition of the geometrician constraints on quantum states and amplitudes, is the crucial point for ensuring the correctness of the model from a discrete geometric point of view, and it is still subject to debate in the literature. In particular, one would expect to find back the Regge action for metric (edge length) variables, or an equivalent classical reformulation, in the semiclassical expansion of the spin foam amplitudes for a generic simplicial complex (or of the corresponding simplicial path integral). Many results are available (see the cited references) on this issue for a single 4-simplex and for extremely simple complexes, at least for the EPRL model, but the results are mixed, and the situation is especially unclear for larger complexes.
The general formula for the spin foam amplitudes, for all the models in this class, in the Riemannian setting, for given cellular complex m dual to the Feynman diagram Γ is the following: The coefficients f are matrix elements of the map between Spin(4) and SU (2) intertwiner spaces: where we have a Spin(4) representation J f labeling each face, a pair of Spin(4) four-valent intertwiners I ve , I v ′ e for every edge and an SO(3) spin j ef for each edge in a given face, while C is the 3j-symbol, and the 15j-symbol is the one of the first type. We have indicated with w the function of group representations that characterizes the implementation of the simplicity constraints defining each model, depending on the representations of Spin(4) and SU(2) labeling each face of the complex. This depends also on the Immirzi parameter γ, through the combination β c−1 c+1 . The amplitude's formula can be rewritten in terms of propagators as follows: where the propagator K is given by: In order to derive the master Integral expression for a GFT Feynman G we have to write down the corresponding (regularized) full amplitude A(G, β, Λ, J → ext ) and then set to zero all the spin associated to the external and contractible internal faces. We are going to give one concrete example of this procedure in the next section, when studying the scaling of the corresponding amplitude.
The models we will deal with in the following are the EPRL model [18] and the Duflo BO model [17], whose defining maps are: where the T function is given by: but it can also be given an expression purely in terms of representation labels. See ref. 17 for more details.
We note that the relative simplicity/complexity of these two models is highly dependent on the basis in which they are expressed, with the flux representation switching such relative complexity with respect to the spin representation given above.
We also point out that other spin foam models, obtained from alternative strategies of imposition of the same geometrician constraints and thus also belonging to the same general class we are considering, can be cast in principle in the same general form, and studied by the same method we will illustrate in the following. Beside difficulties, for some of them, in achieving an explicit and manageable expression for their corresponding w coefficients in representation variables, that makes the analysis more cumbersome, it would indeed be very interesting to perform the same scaling analysis of amplitudes and compare with our results.
More details about the construction of spin foam amplitudes, as well as all the ingredients we mentioned as entering in such construction, in a language well adapted to their GFT embedding, can be found in ref. 17.

RENORMALIZATION OF GROUP FIELD THEORIES FOR 4D QUANTUM GRAVITY
Let us now discuss motivation and current status of renormalization of simplicial GFT (and spin foam) models for quantum gravity.
Beyond the connection to spin foam models and simplicial gravity path integrals, the general strategy for renormalization of GFT models [10, 11,26] is to treat them as ordinary QFTs defined on a Lie group manifold, thus using the group structures (topology, Killing forms, etc) to define 'scales'and mode integration. A natural notion of scale, to be used to label the RG flow, is provided by group representations, which index the spectrum of differential operators on the group, e.g., the Laplace-Beltrami operator, in turn often used to define the propagator of GFT models. Cut-offs imposed as part of a renormalization group scheme are then imposed on representation labels; for example, in the case of SU(2) cutting off the spectrum of the Laplacian operator means imposing the bound d i 1 j i (j i + 1) ≤ Λ 2 , for some real (large) number Λ. This fits well with the fact that divergences in spin foam amplitudes mostly come from the large representations regime. Still, a lot of non-trivial work (beside computational challenges) is needed to adapt for GFTs, whose Feynman diagrams are not graphs but cellular complexes, standard QFT notions, noticing also that any procedure for the contraction of divergent subgraphs of perturbative GFTs has the meaning, from the point of view of the simplicial gravity path integral or spin foam model corresponding to the Feynman amplitudes of the same, of a coarse graining scheme of the corresponding lattice theory.
For a proper renormalization group scheme, however, two more ingredients are needed: control over the theory space corresponding to a given GFT model, i.e. the space of allowed interactions; a detailed characterization of the combinatorics of (the cellular complexes dual to) GFT Feynman diagrams. On neither of these two points much is known for simplicial 4d gravity models. As a result, most work in the context of GFT renormalization has been done focusing on tensorial GFT models, where the above limitations are not present.
Before discussing the goals of GFT renormalization, we spend a few words of caution concerning the physical interpretation of the renormalization group scheme and derived flows. With scales associated to group representation labels, the natural cutoffs entering as UV cutoffs are for large representations. The associated RG then flows from large to small representations (from UV to IR). In LQG and simplicial quantum geometry representation labels identify eigenvalues of geometric operators (e.g., triangle areas or tetrahedral volumes). Large representation labels correspond to large values of such geometric quantities. Thus we have an apparent inversion of roles here, with large distances/volumes playing the role of UV scales in GFT. Caution however should be exercised. In both LQG and simplicial geometry, we know an area of a surface, say, to result roughly speaking from the sum of the individual areas of all elementary surfaces forming the one under consideration, so that one has A ≈ < j > < N > with < j > the average area contribution and 〈N〉 the (average) number of contributions. Moreover, experience from classical Regge calculus and other simplicial gravity formalisms, leads us to expect continuum geometry to be reproduced when the number of elementary excitations contributing to a given continuum geometric quantity is very large, with each contribution smallish (but allowed to be orders of magnitude above Planck size). On the same basis, we expect continuum geometry, and with it any notion of large or small areas, volumes, distances etc, to be the result of coarse graining microscopic, fundamental degrees of freedom like the ones we deal with in the fundamental GFT (or spin foam) formalism. We would better refrain, then, from interpreting simplicial observables directly as geometric, in the sense we attribute to continuum spacetime geometry and physics. Finally, one more alert comes from recalling that, in GFT, the simplicial geometric observables and excitations are the ones associated to the Fock representation of the theory, and probably this is does not correspond to a fully geometric phase in which continuum gravitational (thus, spatiotemporal) physics is to be found, being best adapted to perturbation theory around the fully degenerate (from the point of view of geometry) Fock vacuum.
Also, let us comment on the importance of a better understanding of the symmetry properties of these 4d gravity models, both for the characterization of the corresponding theory space and for their relation to continuum gravity, which is of course crucially characterized by diffeomoprhism symmetry. The issue of symmetries in GFT models is very important but also very much open. At the general level, we do not know much about symmetries of 4d GFT models, beyond the Lorentz invariance of the kinetic and interaction kernels and of the Feynman amplitudes (implemented as in usual lattice gauge theories, since the amplitudes are in fact lattice gauge theories for a (constrained) Lorentz connection). Because of the simplicity constraints and also of the simplicial combinatorics that characterize them, moreover, even the tensorial symmetry typical of tensor models is not present (or at least not manifest). Moreover, even for the few symmetries we know of, in other models, the analysis of their consequences, for example in terms of conservation laws, is complicated by the non-local nature of the GFT interaction (see the analysis [27,28]). Concerning diffeomorphisms, strictly speaking (being smooth transformations) they are not defined in a discrete context like that of GFT Feynman amplitudes, i.e. spin foam models and lattice gravity path integrals, and thus the question becomes whether we can identify some analogue of diffeomorphic symmetry that, in a continuum limit, could be then identified with the one characterizing GR. There are several analyses of such question for 3d (topological) models at the level of spin foam amplitudes [29], lattice gravity (see for example [30][31][32]) and corresponding GFT formulation [33], but nothing similar in the 4d gravity case (where the 4d counterpart of the symmetry identified in the 3d case is actually broken, at the discrete level [34]). When attempting a reconstruction of an effective dynamics of geometry in a continuum approximation, as done in the context of GFT cosmology, one has to proceed in terms of observables of the fundamental theory that have a chance to correspond to diffeomorphic invariant observables in GR, since all the structures of continuum GR on which diffeomorphisms act, e.g., manifold points, directions an coordinate functions, but also fields defined on the same manifold, are simply not present in the theory.

Quantum Consistency and Perturbative Renormalization
GFT models are first defined in perturbative expansion and it is in this perturbative formulation that spin foam amplitudes, and simplicial gravity path integrals, appear. The perturbative GFT amplitudes generically diverge and regularizations have to be imposed. It is this truncation that corresponds to working at a given 'scale'. Is this definition of the quantum dynamics of GFT models consistent? is the spin foam description consistent? Here, consistency means first of all valid for all ranges of dynamical variables, under (controlled) removal of regulators. If not, the GFT model as defined in perturbative expansion, and thus the corresponding spin foam model (and simplicial path integral) cannot be trusted. In the GFT language, this is recognized immediately to be the issue of perturbative renormalizability of a given model. We should only trust, from the spin foam or lattice gravity point of view, only GFT models that turn out to be (perturbatively) renormalizable. We note in passing that there should be no requirement that the model is finite (in the sense of presenting no divergence even before any renormalization); first, we have no obvious reason to expect it, if the model contains an infinite number of degrees of freedom; second, renormalizable models are usually more interesting, as QFTs, than finite ones, since they have a nontrivial RG flow and new effective physics at each scale.
Let us clarify further what we mean, here, to avoid possible misunderstandings. As a general point about field theories, we are saying that finiteness of the Feynman amplitudes associated to a given subset of diagrams, or even to all diagrams involved in a given "scattering process" is not so important, per se, and in fact not necessarily desirable. What is important is that the scattering amplitudes can be -made finiteby suitable renormalization procedure (at any order in perturbation theory), if originally divergent in terms of bare couplings, and after resumming all the diagrams involved in their computation. The final renormalized scattering (or transition) amplitudes are what is physically relevant. A theory that is instead simply finite in the sense of not requiring any renormalization, even if clearly easier to deal with, would be less interesting from a physical point of view because this finiteness would probably indicate that the quantum dynamics is not very rich and it does not change much across scales (i.e. when more of its quantum degrees of freedom are accounted for). The consequence would also be a less interesting phase diagram. This should explain our comments about finiteness of GFT Feynman amplitudes. As for the finiteness of "quantum gravity scattering amplitudes", what we should expect or desire depends on how we interpret the terms. If we take a given GFT model to be a tentative definition of full quantum gravity, then for sure we should hope that its "transition or scattering amplitudes" be finite, in the end, i.e., after renormalization. If the relation between GFT and quantum GR is as in the first case discussed above, i.e., the two are "equivalent", then this also implies that the transition or scattering amplitudes of quantized GR will be finite, when properly defined and after renormalization (even if the renormalization procedure as well as the observables expressing the amplitudes may look very different in the two formulations). If the relation is as in the second case, and thus "quantum GR" is just an effective theory, then we do not have to expect that its transition or scattering amplitudes are finite tout court, but only within the domain of validity of the approximations or truncations leading to it within the fundamental theory.
Thus, the requirement of perturbative renormalizability is an important constraint, which helps removing from consideration inconsistent constructions. Here, the GFT embedding proves potentially very important also for spin foam models (and loop quantum gravity). All known GFT and spin foam models present several ambiguities, some intrinsic to any quantization procedure, others specific to simplicial GFT (and spin foam) models of quantum gravity. Requiring perturbative renormalizability means constraining such ambiguities. To name one, we have little constraints of the face amplitudes of spin foam models, even though they can drastically affect the scaling behavior of the GFT and spin foam amplitudes, to the point of allowing to achieve perturbative finiteness easily by simply fixing them to this end (which also shows why finiteness per se cannot be a goal, without a proper physical understanding), while perturbative renormalizability is a much trickier requirement. For example, see how simple modifications of the Barrett-Crane model (which is also the limit of both EPRL model and Duflo model for infinite Immirzi parameter) affect the resulting amplitudes [35,36].
Let us list some of them. A first one is combinatorial: why restricting to simplicial complexes? These are the ones for which we have a better understanding of the discrete geometry underlying our models, and in particular of the simplicity constraints that characterize them. But what other cellular complexes should be included in the theory for consistency, e.g., because corresponding to the counterterms required for taming the perturbative divergences? Others concern the underlying quantization and imposition of simplicity constraints. Being functions of the flux variables (which are non-commutative), they depend on which quantization map is chosen to quantize such variables. Different choices result in different discrete gravity actions and different simplicial path integral measures, thus different spin foam amplitudes. Also the very definition of the simplicity conditions as operator equations acting on quantum states depends on the chosen quantization map, from which follow thus different constraints on representation variables in the spin foam amplitudes. Further, the strategy by which simplicity constraints are imposed produces in general different models or versions of the same type of models (this is apparent in the Riemannian case, while in the Lorentzian one we only have experience of different versions of the EPRL model). These and other ambiguities are discussed, e.g., in ref. 17. Using SL(2, C) or SU(2) data to label quantum states, which is another choice to make, also leads to potentially different models and amplitudes. Nor one should think that these ambiguities are an artifact of the GFT or spin foam formulation. They can be convincingly argued to be the counterpart of ambiguities in the definition of the canonical Hamiltonian constraint operator and, in a way, failing to fix (at least most of) them via renormalizability conditions would be the counterpart, at the background independent level, of the problem of non-renormalizability of perturbative quantum gravity on a given spacetime geometry [37]. Perturbative GFT renormalizability is thus a crucial issue, also when one looks at it from the perspective of spin foam models, simplicial path integrals or canonical loop quantum gravity.
So, where do we stand, on this important issue? For simplicial GFT models of 4d quantum gravity the answer is, unfortunately, that we are only at the very beginning. The main reasons have been already mentioned. First, we do not know enough of their symmetries to characterize the relevant theory space. Second, the amplitudes for these models are very involved and technically challenging to compute, mostly due to the fact that the imposition of simplicity constraints makes them defined not simply on Lie group manifolds but on particular sub-manifolds of these (usually not even corresponding to homogeneous spaces). Third, dominant configurations (i.e., those giving the most divergent contribution to the amplitudes) are not just flat connections or similarly simple, but correspond to richer configurations from the point of view of simplicial geometry; possibly, they correspond to (or possibly include) the whole set of Regge geometries found as saddle-point configurations in the asymptotic analysis of spin foam amplitudes and corresponding simplicial path integrals. Therefore even power counting results are hard to obtain, and the brute force analysis of divergences is not advanced enough to indicate the needed counterterms, forming the theory space. If the theory space is hard to characterize also in the simpler 3d simplicial case (corresponding to topological BF theory), at least the amplitudes are manageable enough to obtain complete power counting theorems [38], identify some counterterms [39] and nice finiteness results [40].
So computational challenges are one big obstacle. It is on this aspect that we focus in the next section, presenting some new results in the Riemannian context. These new results should be added to other ones we have on the calculation of radiative corrections and basic divergences of both Riemannian and Lorentzian simplicial spin foam models, and on explicit evaluations of their building blocks (mainly the vertex amplitudes). For a partial list, see [41][42][43][44][45][46][47][48][49] and references therein. Future progress will build on these hardwon calculations. In turn these results build on the hard-won understanding, based on both analytical and numerical studies, regarding the asymptotic properties of SU(2) recoupling coefficients and, more recently, of the SL(2, C) recoupling invariants as well, following an extensive ongoing effort in the spin foam community to investigate the behavior of the spin foam transition amplitudes for various models in the constrained BF theory class. A tentative (and incomplete) list of interesting references is refs. 46-refs. 61.
For a comparison, one has to look at the amount of knowledge we have accumulated on tensorial (thus colored) GFT models [10, 11,26]. Here we know several (classes of) models which are rigorously proven to be perturbative renormalizable, comprising both abelian and non-abelian models, on homogeneous spaces, with or without gauge invariance (closure condition), in different dimensions. Divergences are associated to bubbles, i.e. cells of the complex dual to the cellular complex associated to a GFT Feynman diagram, and typically the most divergent diagrams that form the relevant theory space of renormalizable theories are melonic ones, also singled out in tensor models. However, we also know example of TGFT models which are renormalizable outside the melonic truncation [62], and these examples may be relevant also for the case of simplicial GFT models, since the structure of their divergences presents some aspects of the simplicial case.

Continuum Limit and Non-Perturbative Renormalization
GFT models of quantum gravity are bona fide QFTs, thus they possess infinite degrees of freedom, as we expect quantum gravity to do (at least thinking of it naively as a quantum theory of the gravitational field). Control over a very large number of degrees of freedom can only be achieved step by step, within some truncation scheme. With the inclusion of more and more degrees of freedom, we can expect a richer and richer set of new phenomena to be unraveled, simply because the physics of many (quantum, interacting) degrees of freedom is very different from that of few of them. In particular, we expect new phases to be revealed. Controlling the full quantum dynamics is controlling the continuum limit of GFT models, and this implies mapping out as best as we can the phase diagram of the same models. In practical terms, it means being able to evaluate the full GFT partition function, for given values of coupling constants. This is the problem of computing the full non-perturbative renormalization group flow of any given GFT model.
Given the mentioned structural connections, understanding the non-perturbative renormalization of a quantum gravity GFT model implies controlling the continuum limit of the corresponding lattice gravity path integral and spin foam This should already make clear why the precise relation between the (non-perturbative) renormalization of GFT models for 4d quantum gravity and the (non-perturbative) renormalization of continuum quantum GR treated as an ordinary field theory (as in the asymptotic safety approach) can only be envisaged in a very tentative manner. Let us give only some comment on our own tentative perspective on this. take a given GFT model that can be fully defined at the non-perturbative level, thus associated with a continuum phase diagram where RG flow trajectories are welldefined from the deep UV (in the GFT sense) to the full IR (still in the GFT sense), and thus accounting for all the (infinite) degrees of freedom of the model; to achieve this situation is the goal of nonperturbative GFT renormalization, as explained. A matching with GR requires that one can also compute, in the same model, observables which characterize fully a 4d geometry and that can be shown to satisfy the GR equations in a classical approximation. Now, we can envisage two possibilities. If the rewriting is, in the appropriate sense, exact, i.e., if one can in principle go from the GFT formulation of the theory to the geometric "quantum GR" one, in the same continuum limit, then the GFT model could be seen, in fact, as a definition of "quantum GR", without any change in dynamical degrees of freedom. In this case, one could expect that there exist a translation of the RG picture of the given GFT model into the one obtained by a non-perturbative RG treatment of GR, for example as provided (ideally) by the asymptotic safety scenario, and an isomorphism between their corresponding phase diagrams and RG trajectories. If the rewriting requires, instead, some truncation of the dynamical degrees of freedom of the GFT model, is valid only for a subset of the GFT observables, or some other drastic approximation to be valid, i.e., if "quantum GR" turns out to be only an effective, emergent description of some sector of the full quantum GFT, then the situation is different. In this case, we should not expect that the GFT phase diagram matches the GR one, and we can only expect that it will reproduce a portion of it, for scales and regime of couplings where the needed approximations and truncations hold. This regime will probably be the one corresponding to "low energies" from the standard GR and effective QFT perspective. Of course, all the above is very much tentative and it is hard to envisage the precise relation at the current stage of development of GFT as well as of "quantum GR", even though a number of features of GFT models (e.g., the fact that they include a sum over topologies and not just geometries, at least at the discrete level) would suggest that the second scenario is more likely.Where do we stand, at the non-perturbative renormalization level? Beside work on the non-perturbative RG flow of tensor models [67,68], a lot of activity has focused on the analysis of GFT models proper [11,26]. Two main strategies have been followed. One is based on constructive methods, mostly focusing on the resummation of the perturbative series, e.g., showing Borel summability. The other is based on functional renormalization group analysis, either (mostly) based on the Wetterich-Morris equation for the effective action, or the Polchinski equation for n-point functions. For the same reasons that limited work on perturbative GFT renormalization, little is known about the general RG flow of simplicial GFT models of 4d quantum gravity. Simplicial GFT models in 3d have been shown to be Borel summable [69,70] and phase transitions for the GFT formulation of simplicial BF theory in any dimension has been shown to exist [71]. But no similar analysis has been carried over to the 4d gravity case, where, as mentioned, we even lack perturbative indications.
The tensorial GFT case, on the other hand, has been widely explored, mostly via functional renormalization techniques, with many results on a variety of models, again both abelian and nonabelian, with and without gauge symmetries, based on compact as well as non-compact Lie groups, in different dimensions. Concerning UV behavior, asymptotic freedom is found in many examples and asymptotic safety is found in others [72], in various truncations, and the perturbative results have been reproduced from a non-perturbative standpoint. More results on the relevance of truncation beyond the melonic sector have been found [73], and the use of Ward identities for studying the RG flow have been explored [26]. Concerning IR behavior (i.e., the actual continuum limit), work is more limited (and more difficult) at the analytic level, but hints have been found, in various truncations and for various models, of a non-trivial phase diagram. In particular, hints of the existence of Wilson-Fisher fixed points (often found alongside asymptotic freedom in the UV) and of broken (or condensate) phases have been obtained [62,74], indirectly supporting parallel work on the extraction of continuum gravitational physics from such condensate phases [75][76][77].
Even if it is unclear, at this stage, which of these results holds also in the simplicial 4d models, with their additional quantum geometric intricacies, all this work on tensorial GFT models has certainly led to a better understanding of GFT renormalization group schemes (and flows). This will certainly turn out to be useful also for the analysis of full-blown quantum gravity models.

BUBBLE DIVERGENCES AND RADIATIVE CORRECTIONS IN GROUP FIELD THEORIES: SOME NEW RESULTS
In this section we report on some recent results concerning the leading order radiative corrections to N-point functions (N ≤ 6) for the Duflo model and the EPRL model, whose amplitudes we have recalled above. We refer to the cited literature for more details on motivations, construction and features of these GFT (and spin foam) models. Also, we limit our presentation to a summary of results and procedures; a more detailed presentation with be left for a forthcoming publication.

A Warm-Up Example: The 2-Point Function of the Ooguri Group Field Theory Model
As a warm-up, we recall the general procedure to compute the degree of divergences of GFT amplitudes both in the holonomic and spin formulation of GFT models, using the simpler case of the Ooguri GFT model for 4d topological BF theory with local SU(2) invariance (that is, a simplicial GFT model on four copies of SU(2) and with kinetic and interaction kernels made out of delta functions only). Analytical evaluation in group variables -For models defined on the full group manifold or a corresponding homogeneous space (as it is the case in most tensorial GFT renormalization analyses), the evaluation is often conveniently done in the group representation. It can be carried out analytically, and it proceeds as follows. We compute first the bulk or amputated amplitude A bulk by removing all the contributions from the external of the GFT diagram G. This amounts to extracting only the dominant leading divergence of the amplitude (subleading divergences require a more refined procedure). Then we gauge fix the holonomic on all the edges of a maximal rooted tree of the graph G. This reduces the evaluation to involve only a set of gauge-invariant variables. Next we drop the contribution from contractible internal faces 1 (if any). The expression so obtained is the irreducible Master Integral I (G) associated to the amplitude, which can then regularized introducing appropriate cut-offs on the remaining integrals (e.g., by replacing the Dirac-delta functions with heat kernels). It is important to notice that different amplitudes contributing to different correlation functions might reduce to the evaluation of the same master integral. Last we evaluate the remaining integrals. This can be done analytically, exactly or approximately for example via saddle point methods (for example using the abelian asymptotic formula for the heat kernels), to find the Master integral scaling exponent ω(G), i.e. its superficial degree of divergence.
Numerical evaluation in the spin basis -When the analytic evaluation in group variables is not possible, it is often more convenient to pass to the equivalent expression in terms of group representations (like the ones given above for 4d gravity models), and then proceed numerically, along similar steps as in group variables. First we compute the bulk or amputated amplitude by setting to zero all the spins labeling the external faces and using the appropriate identities for degenerate recoupling coefficients. Next we compute the Master Integral I (G), expressed in group representations, by setting to zero the spins labeling internal contractible faces. Its expression can be easily regularized by putting a uniform cutoff Λ on all unbounded summations. Last we numerically evaluate the regularized master integral as a function of the cutoff. This can be done either using the full exact formula or using its approximate asymptotic formula (for large spins) obtained by uniform rescaling of all the spins. The amplitude degree of divergence ω can then be estimated by fitting the data. More precisely it is given by the angular coefficient of the linear best fit in a Log-Log data plot.
Let us illustrate the general procedure with an example. We consider the leading order radiative (melonic) correction to the two-point function of this GFT model. The associated Feynman diagram is shown in Figure 1. The diagram has four external and six internal faces, none of which is contractible.
In the holonomic formulation, the amplitude can be written as follows: In this case, with or without gauge fixing, we can perform all the integrations exactly. Neglecting the contributions from the external faces, without a regularization we would find the divergent result: The master integral can be regularized either via a sharp cut-off or by heat kernels: In both cases the amplitude's degree of divergence reads: The same result can be recovered by evaluating the amplitude in the spin basis. We have: Since there are no contractible faces, by setting to zero all the spins labeling the external faces we immediately obtain the expression of the regularized master integral I Oo . Upon using the appropriate identify for the degenerate 15J-symbol we find: Equation 3.6 is the result of applying the identity A.9 to the Ooguri model's two point amplitude A Oo (G 221 , Λ) (see the combinatorics of the diagram in Figure 1) after setting to zero the spins labeling the external faces of the graph G 221 . Thus the degree of divergence can be obtained by evaluating the master integral's exact formula, or approximately from the above asymptotic formula, by combining the volume factor (replacing the redundant summations) and the face weights with the large-j behavior of the Wigner 6J-symbol, obtaining: ω asy (G melon ) 9 (3.8) 1 An internal face of G is contractible if it has at least one internal edge which is not shared with any other internal face. in agreement with the analytical result obtained in the group formulation.

Radiative Corrections in Simplicial Group Field Theories Models for Quantum Gravity
We now report on some recent results concerning the leading order radiative corrections to N-point functions (N ≤ 6) for the Duflo model and the EPRL model. We identify the relevant 1PI Feynman diagrams, compute the corresponding master integral formulae and use them to evaluate the master integrals' scaling (i.e., the diagrams' superficial degree of divergence) as a function of the cutoff. We also comment on the diagrams' combinatorial properties and on the structure of the corresponding counterterms. Finally we show how these results can be applied 'beyond perturbation theory'to characterize to all orders the scaling of the necklace graphs (an important subclass of diagrams appearing in the radiative correction, also identified as the relevant graphs for renormalizability in the tensorial model of ref. 62). We will derive first general formulae that applied to all models in the chosen class (again, simplicial models constructed from constraining those for topological BF theory) and then specialize to the models of interest by choosing the relevant form for the coefficients w, encoding the geometrician conditions characterizing them.

Leading Order Corrections to the N-Point Functions
The relevant 1PI GFT diagrams, appearing in the perturbative expansion of the 2-point and 4-point functions at the leading  Table 1 with the notation explained in the caption. Let us make a few remarks before moving to the analysis of specific diagrams.
Selected diagrams -The diagrams showed in Table 1 are the only potentially divergent 1PI GFT diagrams at the leading order. All the other leading order corrections to the 3−, 5− and 6−point functions, labeled G 311 , G 511 and G 621 , are manifestly convergent. The diagram G 511 has no internal faces, thus no potentially divergent summation over representations. The diagrams G 311 and G 621 have only one internal contractible face each. Therefore their corresponding Feynman amplitudes are again finite.
2−point diagrams -The diagram G 221 is melonic and therefore also tracial. It has six internal faces and four external ones and it could be expected to be the most divergent LO contribution to the 2-point function. Its dominant (or leading) divergence 2 , if any, can be subtracted by mass renormalization (as done in ordinary QFT). The associated irreducible master integral will be denoted as I 221 . The diagram G 222 has four external faces and four internal faces; one of them, the tadpole face, is contractible. The diagram G 223 has four internal faces and four external ones. Although the diagram G 222 and G 223 are not isomorphic, their Feynman amplitudes can be reduced to the evaluation of the 4−point function's master integral I 421 .
The 4−point diagram -The diagram G 421 is the first melonic correction to the 4-point function. It has three internal faces (forming a bubble) and eight external ones. Its master integral I 421 , associated to the bubble subgraph, controls the UV scaling of the LO non-melonic 2-point diagrams and of all the necklace diagrams contributing to the N-point functions with N 4, 5, 6.
To summarize: in order to determine the scaling behavior and divergent structure of the leading corrections to the 2− and 4−point functions we only need to study the independent master integrals I 221 and I 421 , whose expression we will give and evaluate in the following.

The 2-Point Function
To derive the expression of the master integral I 221 for the leading (melonic) correction to the 2−point function, as explained, we first write down the full regularized amplitude A(G 221 , β, Λ, J → ext ) and then we set to zero all the spins associated to the external faces. Exploiting the identities (A.1, A.10, A.12) and a number of algebraic simplifications, we obtain: Hence after setting to zero the spins J i (j − i , j + i ) labeling the external faces (in this case J 2 , J 5 , J 6 , J 10 according to the labeling conventions adopted in the paper), the four internal propagators will depend only on three spins. Thus in order to derive the Master Integral expression I 221 we first need to compute the formula for the degenerate propagator with one vanishing triad which in turns requires the identity for a type-A NineJ symbols. Upon using the identity A.1 we obtain the degenerate propagator's formula A.12 which together with the A.9 leads us to the Master integral expression 3.9. The same line of reasoning applies to the derivation of the Master Integral formula for any GFT graph with one and only one external leg for each simplicial vertex (like the melonic 2-point function, upon setting the external spins to zero the internal propagators of these graphs will only depend on three spins).
For completeness, the full formula for the propagator is given by: The formula for the degenerate propagator (with one vanishing triad) can be written as follows: Upon setting j ± i j ± the reduced degenerate propagator reads: The above formulas are completely general. For the euclidean EPRL model they can be further simplified yielding the equations  where the (degenerate) propagator K EPRL is given by Equation A.15. A sample of the results of the numerical evaluation of the relevant master integrals 6 is in Table 2.
The degree of divergence ω(G 221 ) is given by the angular coefficient of the linear best fit of the data plotted in logarithmic scale. The mean value 7 of the scaling exponent ω and its standard deviation for the studied cases are summarized in Table 3.
In the case of the Duflo model, some further subtleties arise in the evaluation, due to the more involved nature of the simplicity or geometrician coefficients. These subtleties require additional care in the numerical evaluation of scaling exponents, which are worth emphasizing here, since they are of more general validity in this class of spin foam amplitudes. The master integral formula (3.15) relies on the use of the asymptotic formula for the 9J-symbol (A.8). This might not be very accurate for relatively small spins, like the ones we can concretely explore in our numerical evaluations. Furthermore, since the (A.8) has a stronger suppression rate than other approximate formulas for the 9J-symbol used in scaling analyses, e.g., the equilateral one {9j} ≈ j − 8 3 (in which also all the js are identified, which we cannot do in the Duflo model), the full expression (3.15) is expected to provide only a lower bound for ω. To test this expectation and also to crosscheck the known EPRL results, obtained using the equilateral scaling, we repeat our analysis using in both cases a different asymptotic formula.
In order to derive an asymptotic formula for the melonic master integral I 221 , we localize its expression Equation 3.9 around a background configuration (j − , j + ) by setting j − i j − and j + i j + . For the EPRL model this also implies an homogeneous identification of all the spins j i due to the peculiar form of the EPRL's simplicity coefficients. This procedure has been applied and tested for a number of different simplicial models, including a different version of the EPRL model [41], and it seems to be reliable [44][45][46][47]49]. After appropriate simplifications, the general formula Equation 3.9 becomes.
where μ 1, 2, for the EPRL and Duflo models. The reduced propagator K j − j + is given by Equation (A.14). To determine the amplitude's degree of divergence we combine the scaling of the various factors. The scaling of the equilateral 6J-symbol is given by the Regge formula {6j} ≈ j − 3 2 . According to our analysis, in the large-j regime the propagator (A.14) can be very well approximated by the following expressions: The EPRL formula can be analytically derived, as shown in Supplementary Appendix A (see A. 16, A.17). The corresponding   formula for the Duflo model follows from a direct numerical evaluation of the propagator. It is worth noticing that also the Duflo propagator, containing the much more involved Duflo geometrician coefficients Equation 1.6, peak on the same configurations of Spin(4) representations (j + , j − ). This is due to the asymptotic behavior of such coefficients, which contains also (but not only) the configurations corresponding to the EPRL configurations among the dominant ones [79].
The computed values of the scaling exponent α can be found in Table 4.
The values of ω, obtained by substituting the identities Equation 3.18 into the master integral formula Equation 3.1, are listed in Table 3.
To summarize: the leading order (melonic) correction to the self-energy Equation 3.9 appears to be convergent for the Duflo model and divergent for the EPRL one. The degree of divergence we computed for the EPRL model is in excellent agreement with known analytical results in the literature [41,42]. Concerning the Duflo model for the case l 1 the data clearly indicates that the use of the non-equilateral formula (A.8) in Equation 3.15, as appropriate for this model, strongly suppresses the amplitude scaling leading to a more convergent result. This might also be true for the case l 2 although we cannot state it with full confidence at the moment based on the small cutoff range we tested (recall that Λ max 16 in the full  formula and Λ max 100 in the asymptotic formula). The limited range of Λ values we explored might also explain why the difference between the values of ω for l 1, 2 computed from Equation 3.15 (first and second row of Table 3) is smaller than the difference between the corresponding values obtained from the asymptotic formula (third and fourth row of Table 3).

The 4-Point Function
We now focus on the leading (melonic) correction to the 4point function. The corresponding GFT Feynman diagram G 421 is depicted in Table 1. In order to derive the master's integral expression we follow the same strategy used in the previous section. After setting to zero the spins labeling the external faces and performing the appropriate simplification we find: Once more, the above formula is completely general and valid for any simplicial GFT (spin foam) model for constrained BF theory.For the EPRL and Duflo models it specializes to: where in the first expressions we used the same notations of Equation 3.15.
The degree of divergence of the Master Integral in the Equation 3.20 can be computed by fitting the data shown in Table 5. The resulting values of ω are reported in Table 6.
The scaling of the EPRL 4-point amplitude can be directly read off from the corresponding formula Equation 3. 21.
To summarize: the leading order radiative correction to the 4point function converges for the Duflo model while it diverges quadratically in the EPRL model. Neglecting possible ambiguities in the definition of both models (which, as we emphasized earlier, could affect the face amplitudes and thus the precise scaling behavior) it would then seem that the Duflo model does not require renormalization, at least at this order, while the EPRL model does. But of course higher orders are needed to establish the renormalizability of both models, thus it is hard to draw too many conclusions from this result. More than the divergence degree in itself, it is important to notice that, since the G 421 is melonic (and thus tracial), the corresponding counterterm that is required to absorb the divergence, when present, is proportional to a tensor invariant quartic interaction term (more precisely to the bubble B 41 vertex in Figure 2). Such counterterm is incompatible with a pure simplicial theory space (e.g., a strictly simplicial EPRL model would be non-renormalizable), and this signals the need to extend the theory space of geometric GFT models beyond the simplicial ansatz to include tensorial bubble interactions.

Next-To-Leading Order Corrections to the N-point Functions
We now show how to generalize the master integral formulas for the melonic 2-point function G 221 to an important class of higher order GFT Feynman diagrams.
The master integral expressions Equations 3.9, 3.17 rely on the property that each internal link of the diagram is shared exactly by one external face and three internal ones. Hence, after setting to zero the spins labeling the external faces we are left with: i) a pair of 6J-symbols for each vertex (coming from a pair of degenerate 15J-symbols); ii) propagator of the form (A.13) for each internal edge. A pair of face weights d j − i d j + i for each internal face.
The above combinatorial property is true for any tadpole-free GFT diagram with one and only one external link for each simplicial vertex (here denoted as G NvNvn ). Some examples of these diagrams are shown in Figure 3. Therefore the expression Equation 3.17 can be generalized as follows: with μ 1, 2 respectively for the EPRL and the Duflo model. Upon using the identities Equation 3.18, the asymptotic master integral Equation 3.24 takes the following form: For the Duflo model the value of α must be computed on case by case basis (see Table 4), while for the EPRL model we have α 3 2 . The degree of divergence ω can also be written in terms of the number of vertices N v For the diagrams in Figure 3 the above formula give us the following results, reported in Table 7: the Duflo model amplitudes for all four diagrams are finite; for the EPRL model 8, the first and last diagrams might be logarithmically divergent and therefore require a deeper analysis [78]. The diagram G 442 converges, while the melonic diagram G 441 diverges as a cubic power of the cutoff. The corresponding counterterm is proportional to the quartic bubble vertex B 41 , suggesting again the need for a suitable extension of the theory space to incorporate the relevant tensorial interactions. Once more, this is in fact the crucial lesson we draw from this analysis of divergences, more important, we think, than the precise scaling of the amplitudes, for the reasons already explained.

Beyond Perturbation Theory: The Necklace Diagrams
In the previous two subsections we discussed the all the leading order and a subclass of next-to-leading order radiative corrections to the N-point functions with N ≤ 6. Now we show how to use some of the results we found to estimate the scaling of the so-called necklace diagrams to all orders in perturbation theory.
A (connected) GFT Feynman diagram belongs to the necklace class if (and only if) it consists of an open chain of vertices where each vertex (except the first and last ones) is connected only to its two closest neighbors. Here we restrict to the set of k-necklace diagrams with 4 ≤ k ≤ 6 where k denotes the number of external links, since all other higher order diagrams are either 1-particle reducible or manifestly convergent. Three examples of necklace diagrams are shown in Table 8.
The necklace diagrams share the following remarkable property: the set of internal faces of a k-necklace diagram can always be decomposed into the direct sum of two subsets To summarize: the master integral I 421 encodes the scaling of the necklace diagrams to all orders in perturbation theory. For the EPRL model a consistent (recursive) subtraction of all 8 We point out again that we are studying both models under specific choices fixing the various ambiguities that enter the construction of the spin foam amplitudes. These ambiguities affect, in general, the scaling results. divergences associated to k-necklace graph requires an extension of the theory space to include the appropriate tensor invariant interactions of order four and six. A more detailed analysis of the combinatorial structure of the require counterterms is left for future work [78].

THE ROAD AHEAD: SOME SUGGESTIONS
We close with a set of suggestions for research directions to be pursued, toward a complete understanding the renormalization group flow of simplicial GFT models for 4d quantum gravity.

Scaling and More Scaling, and the Discrete Geometry of Divergent Configurations
The first suggestion is stating the obvious: compute, compute, compute 9 . We need to (pardon us the word pun) scale up the effort in investigating the scaling of simplicial GFT amplitudes, at the same time trying to import insights from tensorial GFT amplitudes. We need to know much more about the divergent configurations and their dependence on the combinatorics of the underlying cellular complex. Lacking better tools, hard bruteforce computations of spin foam amplitudes are the inevitable duty, and these in turn can only build on a better control of the relevant building block the vertex amplitude or, in GFT language, the vertex kernel (which plays an important role also in nonperturbative calculations). Brute force alone will not lead us far, however. On the one hand, we need to develop a more refined analytic understanding of these kernels and resulting amplitudes, and to identify simplified expressions that capture the relevant scaling properties, and their behavior under coarse-graining. On the other hand, where analytical methods do not reach, we need numerical ones to take over; numerical tools for the evaluation of GFT amplitudes are thus badly needed. Packages like have been developed for local QFT Feynmanology would be of course most welcome. On the analytic side, another important objective should be to characterize in detail the (simplicial) geometric meaning of the dominant, most divergent configurations. This is needed to understand the nature of the needed counterterms, but it may also provide insights on the physical features of such GFT models, even beyond their discrete formulation.
The goal here is not to so much to be able to compute GFT Feynman diagrams to arbitrary order (e.g., in vertex or loop expansions). Even in standard QFT, for the physical questions for which the perturbative expansion is the correct approximation scheme, we need to compute (very) many diagrams, but there is often no need nor possibility to go beyond some (usually low) order of approximation and beyond a certain (usually small) number of physical degrees of freedom in the chosen boundary states (of course, the two restrictions go together, since for highly populated boundary states, even the simplest diagrams are of high order). A clear physical picture behind this approximation scheme is as important as computational power. Moreover, coming to the specific quantum gravity case, we would argue that the perturbative GFT expansion, and the description of the dynamics in terms of elementary processes involving few of the fundamental quanta, ie. the usual spin foam language, is not the most convenient approximation to capture the effective continuum physics of quantum gravity, e.g., concerning early cosmology or quantum black holes.
Nor the goal of such analysis of perturbative GFT divergences is establishing that one specific GFT model is finite. Not only renormalizability is a more subtle and possibly interesting feature that finiteness, but one can imagine playing with the ambiguities entering the construction of any given model to modify its scaling behavior and turning it into a finite one. This could be a way to fix or constrain such ambiguities, of course, but it also shows that finiteness per se probably should not be a goal, and that physical conditions fixing the same ambiguities are needed. The main goal of actual computations of GFT amplitudes should then be to provide solid indications on the general power counting of divergences, on the way to a renormalizability proof, and, even more, to indicate the relevant counterterms to be added to the model, and thus the relevant theory space of the starting GFT model. More generally, the goal of such perturbative calculations should be to provide information and tools to be employed to go beyond the perturbative setting and dwell into non-perturbative GFT renormalization. It is only the latter that can provide us with the insights and the results we need to truly explore the formal solidity and effective continuum physics of GFT models of quantum gravity.

Group Field Theories Theory Space, Colors and Relation Between Simplicial and Tensorial Models
We have emphasized several times already the importance of defining the relevant theory space of simplicial GFT models, in order to set up a proper renormalization scheme (perturbative and non-perturbative). Much more work should be devoted to this issue, in particular understanding more about the symmetries of such quantum gravity models. One question is whether the yet to be identified theory space of simplicial GFT models relates to the one of tensorial GFTs. We speculate that they do and, in fact, some hints that the two may largely coincide are known. First, taking seriously the tensorial nature of GFT fields implies coloring (thus distinguishing and ordering) their arguments and, as a consequence, their Feynman diagrams. As 9 To be clear: 1-loop and 2-loop calculations would be probably enough to extract a lot of interesting properties from the perturbative expansion of GFT models for 4d quantum gravity. They may even be sufficient, since we actually expect that most interesting physics should be looked for elsewhere, i.e., within different regimes and approximations of the fundamental quantum dynamics, and not in the perturbative (spin foam) expansion. Moreover, as we stressed, such perturbative calculations should have as main goal to identify the theory space within which the relevant GFT models should be placed, more than simply pushing the perturbative analysis for its own sake. We emphasize, however, that, in fact, a lot of computational effort is required to perform such 1-loop and 2-loop calculations because one needs to consider many and complicated GFT diagrams, and involved spin foam amplitudes, already at this order. These are the calculations we intend to encourage. we noted, this coloring allows a precise control over the topology of the cellular complexes dual to these Feynman diagrams [7] and, in turn, this greater control allowed for many results that are central in renormalization analyses of tensorial GFTs (e.g., large-N expansions) [10,11]. A precise control over the topology of the Feynman diagrams, i.e., the cellular complexes on which spin foam amplitudes are based, is needed also in simplicial GFTs, if one aims at identifying the nature of divergences, leading to precise power counting results. It is also needed for identifying key symmetries, as we know already in the case of topological BF models [80], where the complete power counting also relied on the full topological information on the underlying cellular complex [38]. Thus we have a strong argument for relying on colors also in simplicial GFT models of 4d quantum gravity; the form of the corresponding spin foam amplitudes would remain unchanged, but they would now be defined on full 4d cellular complexes, rather than just their 2skeleton. Assuming we work on colored simplicial GFT models, we then have two preliminary results that suggest a close relation with tensorial models. One is that using colors one can identify similar symmetries in the simplicial case than one finds in the tensorial one [28]. The other is that integrating out all fields except one in a colored simplicial GFT model (in any dimension, with trivial kinetic term) produces an equivalent tensorial GFT model for the remaining field (with the same coupling constant for all interactions) [81]. A third general fact pointing in the same direction is that divergences in simplicial GFT models for topological BF theory, which is the starting point of the construction of simplicial 4d gravity models, are associated to bubbles in the cellular complex, which are in fact the cells associated to allowed interactions in tensorial GFTs with the same base group manifold. These results, in our opinion, suggest that there could be a single theory space containing both (colored) simplicial models and tensorial ones, with interaction kernels in the tensorial directions yet to be identified.

Group Field Theories Models With Local Directions
The third suggestion for further research is to devote attention to the renormalization of GFT models which combine the combinatorially non-local pairing structure on geometric variables, in GFT interactions, with the presence of local directions. This includes both simplicial GFT models and tensorial ones, with the distinction referring to the pairing of geometric variables.
There are two main examples of such 'mixed'models. One is the tensorial models used to describe SYK-like many-body systems [82], whose renormalization has been in fact studied in several cases. Here the non-local, tensorial indices are usually reduced to finite sets (we have thus simple tensor models, rather than full GFTs) and the single local direction is a time variable. The standard SYK models are indeed quantum mechanical models in 0 + 1 dimensions, with generalizations to higher dimensions (thus, with more local directions) having been proposed. The other class of mixed models is the extension of (simplicial) GFT quantum gravity models to include scalar fields coupled to gravitational degrees of freedom [83]. These extended models have been studied in particular in the context of GFT condensate cosmology [75-77, 84, 85], with the additional scalar fields playing (also) the role of clock and rods that allow to define relational, diffeo-invariant observables in terms of which an effective cosmological dynamics can be extracted from the GFT hydrodynamics.
The potential physical interest of these models, and of their renormalization analysis, is thus obvious. They present several interesting issues. The presence of both local and non-local directions may modify sensibly the renormalization flow and the structure of divergences, thus leading to different dominant diagrams and effective dynamics in both UV and IR sectors. One can also envisage setting up an altogether different renormalization group scheme, adopting a notion of scale tied to the scalar (local) directions, rather than the group manifold (or involving both), potentially producing very different results. Such focus on the flow parametrized by variables with a (tentative) physical interpretation as relational time/space variables may also allow a more direct physical interpretation of the renormalization flow itself, e.g., in a cosmological context (even though similar cautionary remarks as for the usual renormalization scheme would apply here).

Relation With Lattice Spin Foam Renormalization
We have emphasized how renormalizing a GFT model is tantamount to renormalizing (and studying the continuum limit of) the corresponding discrete gravity path integral and spin foam amplitudes, from a different standpoint. But the GFT formalism is only one way to provide a complete definition of spin foam models, the other being to view them as a peculiar (because background independent) lattice theory and setting up some appropriate refinement procedure. Therefore, it would be very important to compare results obtained in the context of GFT renormalization, especially for simplicial quantum gravity models, with the results and techniques developed for renormalizing spin foam amplitudes from a lattice gauge theory perspective [65,66,[86][87][88].
In this lattice-focused approach to spin foam renormalization, a cut-off is also imposed on representation variables, but the notion of 'scale'is rather given by the combinatorial complexity of the underlying lattice, and the renormalization group flow is driven by refinement/coarse-graining steps ordered by such complexity. Refinement/coarse-graining steps affect both bulk lattices and boundary graphs, and the flow of quantum amplitudes is constrained by the requirement of their consistency under restriction to coarser boundary states.
Despite their differences, the two renormalization schemes share several, since also GFT subtraction moves amount to lattice coarse-graining steps, and corresponding maps between associated amplitudes are also built-in in the (perturbative) QFT renormalization steps used in the GFT context. Still, a detailed work of translation between the two frameworks would be very useful. This work may require, on the GFT side, a combination of functional renormalization group techniques, since we are interested in the continuum limit of spin foam models, and perturbative expansions, given that spin foam models arise in such expansion. This comparison would be beneficial for both approaches; in particular, it would emphasize the role of combinatorial complexity of boundary states in the GFT renormalization flow. This work should be carried out for all models that have been studied in both settings (also in the lattice renormalization approach work has been confined mostly to highly simplified models), aiming of course at unraveling the continuum phase diagram of 4d quantum gravity from two perspectives at once.

Group Field Theories Renormalization Via Tensor Networks
One powerful set of techniques coming from the theory of quantum many-body systems, that have been already applied in the context of lattice-based renormalization of spin foam models, uses the language of tensor networks [89,90]. This language is useful both for numerical studies and for emphasizing the role of entanglement in the renormalization group flow [91,92]; in particular, it allows to unravel topological quantum phases of many-body systems.
In the case of GFT models, the interest in importing techniques from tensor networks goes beyond these general facts, and stems also from the fact that GFT states themselves can be seen as generalized tensor networks [93], and by the related fact that entanglement is responsible for the basic connectivity between GFT quanta that gives rise to extended discrete structures labeled by quantum geometric data. The many facets of the GFT formalism, moreover, would allow for a manifold application of tensor network techniques. On the one hand the basic GFT field is a tensor and its quantum states are tensor networks, as mentioned; on the other hand, it remains a QFT, calling for continuum tensor network techniques as employed, say, in standard scalar quantum field theory [94]. At the same time, its Feynman amplitudes are lattice gauge theories, to which a different set of tensor network techniques can be applied [95] (as developed in the context of spin foam lattice renormalization). And they remain quantum many-body systems, peculiar for their background independent nature, but still conventional enough to allow the deployment of tensor network methods taken from their natural context.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.