Pietro Dona
email: dona@cpt.univ-mrs.fr

Francesco Gozzini
email: gozzini@cpt.univ-mrs.fr

Giorgio Sarno
email: sarno@cpt.univ-mrs.fr

Pietro Donà

Numerical analysis of spin foam dynamics and the flatness problem

published or not. The documents may come

Numerical analysis of spin foam dynamics and the flatness problem

Introduction and Motivations

The EPRL-FK spin foam theory [START_REF] Engle | LQG vertex with finite immirzi parameter[END_REF][START_REF] Freidel | A New Spin Foam Model for 4D Gravity[END_REF] is an attempt to define the dynamics of loop quantum gravity. It provides a regularized, background-independent and Lorentz covariant quantum gravity path integral with the definition of a partition function on a triangulation. The theory assigns transition amplitudes to spin network states on the boundary of a triangulation.

Recovering General Relativity (or at least a discrete version of it) in the semiclassical limit is a fundamental test for spin foam theories. The large spin limit of the EPRL vertex amplitude has been largely explored in both euclidean [START_REF] Barrett | Asymptotic analysis of the eprl four-simplex amplitude[END_REF] and lorentzian versions [START_REF] Barrett | Lorentzian spin foam amplitudes: graphical calculus and asymptotics[END_REF]. Remarkably, under a homogeneous rescaling of all the spins, the single vertex amplitude with coherent boundary data contains the Regge action, a discrete formulation of General Relativity.

However, in the case of extended triangulations, the semiclassical limit of the theory is not completely understood. In particular, the question of what semiclassical geometries dominate the summations over bulk degrees of freedom is still open.

Different calculations suggest that the summation over bulk degrees of freedom is dominated by flat geometries 1 . This observation, often called flatness problem, has been mentioned by Freidel and Conrady [START_REF] Conrady | On the semiclassical limit of 4D spin foam models[END_REF], Bonzom [START_REF] Bonzom | Spin foam models for quantum gravity from lattice path integrals[END_REF], and Han [START_REF] Han | Semiclassical Analysis of Spinfoam Model with a Small Barbero-Immirzi Parameter[END_REF]. Hellmann and Kamiński [START_REF] Hellmann | Holonomy spin foam models: Asymptotic geometry of the partition function[END_REF] proposed a different analysis, based on the calculation of the wavefront set of the spin foam partition function and a geometrical interpretation of its variables, and they found a similar problem.

A more refined stationary phase analysis of the euclidean EPRL model has been carried out by Oliveira in [START_REF] Oliveira | EPRL/FK asymptotics and the flatness problem[END_REF]. They include all the constraints on the boundary data necessary to obtain a non-exponentially suppressed vertex amplitude. The initial claim that the flatness problem is not present in that setting has been recently questioned by Kaminski and Engle [START_REF] Engle | private communications and ILQGS seminar[END_REF]. It has been suggested that taking only large spins is not the right way to access the semiclassical regime of the theory [START_REF] Magliaro | Emergence of gravity from spinfoams[END_REF][START_REF] Han | On Spinfoam Models in Large Spin Regime[END_REF][START_REF] Asante | Effective Spin Foam Models for Four-Dimensional Quantum Gravity[END_REF].

The importance of numerical methods for the investigation of spin foam theories grew considerably in recent years. The software library sl2cfoam is a C based high-performance library and is the foundation for the numerical computation of general spinfoam amplitudes, both in 3 and 4 dimensions, with topological BF or lorentzian EPRL models [START_REF] Dona | Numerical methods for eprl spin foam transitionamplitudes and lorentzian recoupling theory[END_REF]. The library is based on the factorization [START_REF] Speziale | Boosting Wigner's nj-symbols[END_REF] and was used to explore the large spin limit of both models [START_REF] Dona | SU(2) graph invariants, regge actions and polytopes[END_REF][START_REF] Dona | Numerical study of the lorentzian engle-pereira-rovellilivine spin foam amplitude[END_REF].

While still at the early stages, the evaluation of many-vertices spinfoam amplitudes is possible. In a recent work [START_REF] Dona | Searching for classical geometries in spin foam amplitudes: a numerical method[END_REF], we have developed an algorithm to determine the existence and estimate the position of stationary phase points in the summations over the spins of internal faces. We tested its effectiveness applying it to BF spin foam theory in 3D (the Ponzano-Regge model), where the stationary phase points are directly connected to the solution of the equations of motion of euclidean Regge calculus. We adapt the algorithm to the BF spin foam theory in 4D. This topological model is the starting point for the construction of physical spin foam models like the EPRL model. For this reason, there are many similarities between the two. For example, the boundary states of the transition amplitudes are the same, as well as the semiclassical geometries emerging from the single vertex asymptotics. What geometries dominate the BF transition amplitudes in the large spin limit is an interesting question on its own.

In this work, we focus on the ∆ 3 triangulation, formed by three 4-simplices sharing a common triangle. It is the simplest triangulation with many vertices and a bulk face, and it is the standard example used in the flatness problem literature [START_REF] Bayle | Short-scale Emergence of Classical Geometry in Euclidean Loop Quantum Gravity[END_REF][START_REF] Oliveira | EPRL/FK asymptotics and the flatness problem[END_REF][START_REF] Hellmann | Holonomy spin foam models: Asymptotic geometry of the partition function[END_REF].

Solutions of the classical equation of motion dominate the summation over the bulk degrees of freedom in the large spin limit. We expect them to manifest through the presence of stationary phase points. Note that in classical Regge calculus, the ∆ 3 triangulation is too simple to give non-trivial dynamics since boundary conditions fix all the lengths. Hence, this triangulation is not a good example to probe the large spin limit of the EPRL model if we suppose the classical underlying theory to be (area-angle) Regge Calculus. Nevertheless, the study of this amplitude can teach us some valuable lessons.

We find an unexpected result. We construct three examples of the 4D euclidean Regge ∆ 3 triangulation. One is the flat triangulation characterized by a zero deficit angle around the shared triangle. The other two are curved triangulations. We compute the coherent amplitude associated with each triangulation, and we look for stationary phase points in the internal spin, dual to the shared triangle. We find the presence of a stationary phase point in all the cases. We estimate its value, and we find it compatible with the area of the dual triangle of the prescribed triangulation.

If we apply arguments similar to the ones claiming flatness of the EPRL model to BF theory, we would infer that curved triangulations should be suppressed in the sum. Therefore we would not expect to find any stationary phase point in this case. This shows tension between our results and the arguments declaring that the EPRL model is flat.

The paper is organized as follows. In Section 2 we review the spin foam formulation of BF theory and the ∆ 3 transition amplitude. Section 3 describes the Regge triangulations we use to construct the boundary data of the amplitude. In Section 4 we illustrate the algorithm we use to search for stationary phase points and its application. We review the main arguments of the flatness problem of the EPRL model applied to the topological BF model in Section 5. We conclude with a summary of our numerical results and their implications.

All the computations are performed on the CPT servers in Marseille. Each machine is equipped with a 32 cores CPU Intel(R) Xeon(R) Gold 6130 with a base frequency of 2.10 GHz and 196 GB of RAM. In the public repository [20] we publish the C code used to calculate the transition amplitudes along with the Wolfram Mathematica notebooks used to prepare the boundary data and to perform the numerical stationary phase point analysis.

BF theory and the ∆ 3 transition amplitude

General Relativity in four dimensions can be written as a BF theory with constraints. This is the starting point in the construction of many spin foam models, as the lorentzian EPRL model. In this paper, we study the spin foam formulation of the SU (2) BF theory on a four-dimensional manifold. For a detailed overview of BF theory and its relations with spin foam models, see [START_REF] Baez | An Introduction to spin foam models of quantum gravity and BF theory[END_REF] or [START_REF] Perez | The spin foam approach to quantum gravity[END_REF]. A triangulation ∆ of the manifold is dual to a two-complex ∆ that consists of a set of vertices (dual to the 4-simplices), edges (dual to the tetrahedra) and faces (dual to the triangles). We assign an SU [START_REF] Freidel | A New Spin Foam Model for 4D Gravity[END_REF] holonomy g e to each half-edge of the triangulation. The partition function of the theory discretized on ∆ is given by

Z(∆) = e dg e f δ(g e1 . . . g en ) , with e i ⊂ f , (1)
the product g e1 . . . g en is the holonomy around the face f , and dg e is the Haar measure on SU (2). We expand the delta function using the Peter-Weyl theorem as δ(g) = j (2j + 1)T r(D j (g)), where D j (g) are the Wigner matrices of the SU(2) representation of spin j. We get

Z(∆) = j f e dg e f (2j f + 1)T r(D j f (g e1 . . . g en )) (2)
=

j f f (2j f + 1) v . (3)
It is useful to introduce a graphical notation for spin foam vertices. Every box represents one integral over SU [START_REF] Freidel | A New Spin Foam Model for 4D Gravity[END_REF], each line is an irreducible representation of SU(2), and each edge contains four lines. For more details on graphical notations see [START_REF] Varshalovich | Quantum Theory of Angular Momentum[END_REF][START_REF] Dona | SU(2) graph invariants, regge actions and polytopes[END_REF]. We reorganized the partition function as a product over the vertices of ∆ (4-simplices of ∆) and a sum over all the spins. Furthermore, the integrations over SU(2) can be performed explicitly, introducing a sum per edge over a 4-valent intertwiner space i e . We obtain the familiar form for the partition function in terms of amplitudes

Z(∆) = j f ,ie f A f e A e v A v . (4)
The face amplitude and the edge amplitude are given respectively by the dimensional factors A f = 2j f + 1 and A e = 2i e + 1. The vertex amplitude A v is given by a {15j} symbol. In this work we use the {15j} symbols of the first type (see Appendix A for our conventions). Curiously, the vertex amplitude of this theory with coherent boundary data contains the Regge action of a 4-simplex in the large spin limit. This connection with discrete gravity has been extensively studied both analytically [START_REF] Barrett | Quantum gravity asymptotics from the su(2) 15j symbol[END_REF] and numerically [START_REF] Dona | SU(2) graph invariants, regge actions and polytopes[END_REF], and is limited to a single vertex. The topological nature of the model is encoded in the delta functions that appear in triangulations with internal faces.

We study the spin foam transition amplitude of the ∆ 3 triangulation. The 2-complex dual to the ∆ 3 triangulation is represented in Figure 1. It is the simplest amplitude with 3 vertices, 18 boundary faces and one internal face. Each vertex shares an edge with the other two and has three boundary edges.

The boundary data consist of coherent states parametrized by 18 spins and 36 unit normals, one spin and two normals per each strand of the spin foam diagram. We denote with j abc each boundary spin and with x the internal spin. In the next section, we study the four-dimensional geometry we use to specify the boundary data; the notation will be clearer then. The spin foam transition amplitude, omitting the summation over x and the associated face amplitude, is represented using the graphical notation in Figure 1. We refer to [START_REF] Hörmander | The Analysis of Linear Partial Differential Operators I: Distribution Theory and Fourier Analysis[END_REF] in appendix C for the complete picture with all the spin labels and orientations of the faces.

x We perform the SU(2) integrals and obtain the ∆ 3 spin foam amplitude in terms of {15j} symbols (see appendix C for more details)

W ∆3 (j f , n f a ) = (-1) χ x w ∆3 (x, j f , n f a ) (5)
= (-1) χ x (-1) x (2x + 1) We associate to each boundary tetrahedron a coherent state characterized by 4 spins, the areas of its faces, and 4 unitary vectors, the outward normals to the faces of the tetrahedron. The triangulation is made of points, segments, faces and tetrahedra. We use the triangulation representation in Figure 2 as reference. The boundary spins j abc are labeled by a triple of points a, b, c since they correspond to triangles. We decide to focus on Regge geometries, fully characterized by their 15 lengths ab of the segment joining points a, b. Notice that all the lengths belong to at least one boundary tetrahedron, therefore fixing the boundary determines the triangulation completely.

k b b (-1) k b (2k b + 1) ia a (2i a + 1)c ia ( n f a ) ×    i 1 j
For simplicity and easier computability, we restrict our analysis to geometries with a high degree of symmetry. We require "cylindrical" symmetry forcing the three 4-simplices to be identical. This symmetry requires the bulk triangle to be equilateral ( 13 = 35 = 15 ) and the bulk tetrahedra to be isosceles, having an equilateral base and three identical isosceles faces ( 12 = 25 = 23 = 16 = 36 = 56 = 34 = 14 = 45 ). The remaining three lengths also have to be equal ( 26 = 24 = 46 ). Furthermore, we require all the areas of the boundary triangles to be the same, and this fixes 13 = 26 . The symmetry reduces the degrees of freedom of the Regge triangulation from 15 to 2, for example, 12 and 13 . In terms of the lengths, we can compute all the geometric quantities of the triangulations: areas, three and four-dimensional volumes and two, three and four-dimensional dihedral angles. It is interesting to compute the area of the boundary triangles λ and α the 4D dihedral angle around the bulk triangle of one 4-simplex:

λ = 1 4 13 4 2 12 -2 13 , sin 2 α = 3 4 2 13 12 2 12 -7 2 13 ( 2 13 -3 2 12 ) 2 . (6)
They are independent variables and under a linear rescaling of the lengths λ scales quadratically while α does not scale. In the following, we identify with λ the scale of the triangulation. Another interesting geometrical quantity for our analysis is the value of the area of the shared triangle

x g = √ 3 4 2 13 = 3 6 + sin 2 α + 6 1 -sin 2 α 48 + sin 2 α 1 2 λ . (7)
For Regge geometries curvature can be expressed in terms of deficit angles associated to triangles shared between 4-simplices. In the particular class of symmetric triangulations we are considering the deficit angle associated with the shared triangle is

δ = 2π -3α . (8)
We parametrize the triangulations in this paper with λ and δ instead of the two lengths. Once the triangulation is specified, we compute all the normals of the boundary tetrahedra to use them in defining the boundary coherent states.

In this work we consider the flat triangulation characterized by δ = 0 at various scales λ, a first curved triangulation characterized by δ ≈ 3.60 at scale λ = 20, and a second curved triangulation characterized by δ ≈ 2.47 at scale λ = 20. The chosen parameters of the curved triangulations are particularly convenient for this numerical analysis. We motivate our choice of focusing on a single scale for the curved configurations in the next section. All the notebooks we used to compute all the geometric data and the boundary normals for the three triangulations are available in [20].

Numerical analysis

The algorithm in [START_REF] Dona | Searching for classical geometries in spin foam amplitudes: a numerical method[END_REF] is designed to determine the presence and estimate the value of stationary phase points in summations over spin foam bulk spins. We adapt it to the amplitude in examination that we study for different choices of boundary data. In this section, we briefly recall features of the algorithm, and we highlight the differences with the original version. We refer to the original paper for a detailed explanation.

The amplitude (5) involves a sum over the spin x of the internal face and the intertwiners k b of the three internal edges. For brevity let us omit the dependence on the boundary data j f , n f a on the summand w ∆3 (x) ≡ w ∆3 (x, j f , n f a ). First, we compute w ∆3 (x) for each value of the internal spin x using the library sl2cfoam [START_REF] Dona | Numerical methods for eprl spin foam transitionamplitudes and lorentzian recoupling theory[END_REF]. We focus on two quantities

P w (x) = x x = xmin w ∆3 (x ) , and R w (x) = x ∈I c x w ∆3 (x ) , (9)
where I c

x is an interval centered in x with width 2c. We call them respectively the partial sum and the running sum. The algorithm is based on the behavior of these two quantities in the presence of stationary phase points for W ∆3 .

Suppose that x 0 is a stationary phase point for W ∆3 . The sum over the internal spin interferes destructively away from x 0 and constructively near x 0 . The partial sum P w (x) stays roughly constant away from x 0 and has a sudden "jump" near x 0 . Similarly, the running sum R w (x) is close to zero away from x 0 and peaks at x 0 . The width of the peak depends on the parameter c, the size of the interval that characterizes the running sum.

Second, we consider R w (x) = |R w (x) c1 R w (x) c2 R w (x) c3 | correlating three running sums with different interval sizes c i and we use a Mathematica's built-in function to locate the peaks of R w (x). We repeat this step for all the possible triples of c i in a reference set. Since we are working at smaller spins then [START_REF] Dona | Searching for classical geometries in spin foam amplitudes: a numerical method[END_REF] in this work we use as reference set c = √ x max -x min /2 plus or minus 50%. The numerical estimate of the stationary phase point is given by the average of the outcomes of this procedure with an error given by the standard deviation (both values are rounded to the nearest integer). This aims to eliminate the dependence of our analysis from the choice of a particular interval.

The first triangulation we study is the flat one. We report in Figure 3 the running sum and the partial sum for this configuration at λ = 30.

••••••••••••••••••••••••• • • • • • • • • • • • • • ••••••••••••••••••••••••••••••••••• ••••••••••••••••••••••••••••• • •• •• • •• • • • • • • • •• • • ••••••••••••••••••••••••••••••••••••••••••• Figure 3: Flat configuration.
Left: Averaged running sum Rw(x) with interval sizes 2,6, and 10. We normalize respect to the value of its maximum. Right: Real part of the partial sum Pw(x) normalized respect the absolute value of the amplitude. In both plots, we mark with a red line the value of the geometrical area of the bulk triangle, and with a green band the values of the numerical estimate of the stationary phase point within one standard deviation.

Our algorithm is implemented in Mathematica, the notebooks are available in [20]. We estimate the presence of a stationary phase point for spin x n0 = 40 ± 1. The value of the geometric area shared by the three 4-simplices [START_REF] Han | Semiclassical Analysis of Spinfoam Model with a Small Barbero-Immirzi Parameter[END_REF] for the flat configuration is x g0 = 1.34 λ ≈ 40.2. It is in perfect agreement with our numerical analysis.

To corroborate the result we repeat the analysis at other scales λ = 10, 12, 14, 16, 18, 20. We summarize in Table 1 the numerical estimate of the stationary phase point and the geometric area of the shared triangle. Signatures of the stationary phase point are present at all scales, and its position is always compatible with its geometric counterpart. We also note the resources necessary to perform the computations on our server. numerical analytical memory (GB) time (h) λ = 10 :

13 For each calculation, we report the time (in hours) and the memory (in GB) required for these computations on our server. Time and memory for λ = 30 are an estimation based on the previous points. The computation, in this case, was performed on multiple machines, a piece at a time and took more than two months to complete.

±
Notice that the agreement is excellent even at small scales. This confirms the result of [START_REF] Dona | SU(2) graph invariants, regge actions and polytopes[END_REF][START_REF] Dona | Numerical study of the lorentzian engle-pereira-rovellilivine spin foam amplitude[END_REF], the asymptotic regime is reached already at relatively low spins. Moreover, at larger scales λ, the computation becomes computationally demanding both on memory and time. This motivated us to focus only on λ = 20 for other boundary data, the best compromise in terms of results and our computational resources.

We perform our numerical analysis also for boundary data corresponding to the two curved triangulations we described in 3 at scale λ = 20. We report in Figure 4 and 5 the running sum and the partial sum for these configurations.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •••••••••••••••••••• •• • •• • •• • • • •• • •••••••••••••••••••••••••••
Figure 4: First curved configuration. Left: Averaged running sum Rw(x) with interval sizes 2,6, and 10. We normalize respect to the value of its maximum. Right: Real part of the partial sum Pw(x) normalized respect the absolute value of the amplitude. In both plots, we mark with a red line the value of the geometrical area of the bulk triangle, and with a green band the values of the numerical estimate of the stationary phase point within one standard deviation. In both cases, we see a clear signature of the existence of a stationary phase point. For the first curved configuration we estimate the position of the stationary phase point for spin x n1 = 28 ± 1. It is compatible with the geometric area (7) of the curved triangulation x g1 = 1.39 λ ≈ 27.8. Analogously, for the second curved configuration we find a stationary phase point at spin x n2 = 25 ± 1. Also in this case it is compatible with the expected geometric area x g2 = 1.26 λ ≈ 25.2.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •••••••••••••••• • • • •• •• • • • •• • •• • •• •••••••••••••••••••••••••••

Flatness arguments applied to W ∆ 3

In this section, we repeat the principal arguments for flatness of the EPRL model, on BF theory. To be specific we focus on the W ∆3 amplitude [START_REF] Conrady | On the semiclassical limit of 4D spin foam models[END_REF]. We sketch the various arguments to keep the exposition as clear as possible, and we refer to the original papers for more details. The arguments [START_REF] Conrady | On the semiclassical limit of 4D spin foam models[END_REF][START_REF] Bonzom | Spin foam models for quantum gravity from lattice path integrals[END_REF][START_REF] Han | Semiclassical Analysis of Spinfoam Model with a Small Barbero-Immirzi Parameter[END_REF] are based on the large spins approximation of the vertex amplitude. They approximate the integrals over the group using saddle point techniques, and then they look for stationary phase points in the sum over the bulk degrees of freedom. Inserting a resolution of the identity in terms of SU(2) coherent states in each internal edge the summand of the amplitude reduces to

w ∆3 (x, j f , n f a ) ∝ A 1 (x, j f , n f a )A 2 (x, j f , n f a )A 3 (x, j f , n f a ) (10)
where A i (x, j f , n f a ) are the {15j} symbols contracted with coherent intertwiners. The proportional symbol in [START_REF] Engle | private communications and ILQGS seminar[END_REF] indicates that we are omitting multiple dimensional factors and the integrals coming from the resolution of the identities we inserted. They both play no role in the calculation we illustrate in this section.

If the boundary spins j f are large, we can safely assume that also the spin of the bulk face x is large. The asymptotic expansion of the coherent {15j} symbols was studied in [START_REF] Barrett | Quantum gravity asymptotics from the su(2) 15j symbol[END_REF][START_REF] Dona | SU(2) graph invariants, regge actions and polytopes[END_REF]: if the boundary data of the coherent vertex forms the boundary a euclidean 4-simplex, the amplitude can be approximated with

A i (x, j f , n f a ) ≈ N i (x, j f , n f a ) cos(S i (x, j f , n f a )) (11)
where the factor N i contains numerical constants and the hessian, both of them are not relevant for the following considerations. The function S i = f ∈vi j f θ f,i + xθ x,i is the Regge action for the 4-simplex: the sum runs over the triangles f (faces) belonging to the ith 4-simplex (vertex) v i , j f are the areas of the triangles and θ f,i the 4D dihedral angles around the triangle f of the 4-simplex i. The dihedral angles can be reconstructed from orientation invariant scalar products between the normals n f a using the spherical cosine law. In this analysis, we assume that vector geometries play no role or can be ignored by selecting the boundary data for the ∆ 3 amplitude appropriately.

In the large spin regime ( 10) is proportional to

w ∆3 (x) ∝ cos(S 1 (x)) cos(S 2 (x)) cos(S 3 (x)) = e i(S1(x)+S2(x)+S3(x)) + e i(S1(x)+S2(x)-S3(x)) + ( 12
)
e i(S1(x)-S2(x)+S3(x)) + e i(-S1(x)+S2(x)+S3(x)) + c.c. ,

where we omitted the dependence on boundary spins j f and normals n f a that we consider fixed from this point on. By linearity, we can search for stationary phase points of each of the eight terms in [START_REF] Han | On Spinfoam Models in Large Spin Regime[END_REF] independently and sum the results. We focus on the first term, the analysis of the others is similar. The derivative of the Regge action S i respect to the area x is given by δS i (x)/δx = θ x,i . The stationary phase points are the solutions of the equation

δ δx (S 1 (x) + S 2 (x) + S 3 (x)) = θ x,1 + θ x,2 + θ x,3 = 0 (13)
and this apparently implies that flat geometries dominate the summation over the bulk degrees of freedoms in the large spin limit. Notice that if we complete the analysis with the other terms in [START_REF] Han | On Spinfoam Models in Large Spin Regime[END_REF] we would obtain that the oriented deficit angle is vanishing. In this section we disregard this additional but distinct complication, that was discussed in length in the proper vertex literature [START_REF] Engle | Proposed proper Engle-Pereira-Rovelli-Livine vertex amplitude[END_REF].

The main counterargument to this analysis is that (11) holds only if some constraints (closure and shape matching) are satisfied by the boundary data. Therefore, the variables are not independent and we cannot take the variation in [START_REF] Asante | Effective Spin Foam Models for Four-Dimensional Quantum Gravity[END_REF] without including the constraints in the action.

Recently Oliveira [START_REF] Oliveira | EPRL/FK asymptotics and the flatness problem[END_REF] proposed a more accurate stationary phase point analysis taking into account also these constraints. They studied the spin foam transition amplitude on the ∆ 3 2-complex in the euclidean EPRL model. The original paper claimed that in this setting the flatness problem was not present, however Kamiński and Engle [START_REF] Engle | private communications and ILQGS seminar[END_REF] corrected an oversight obtaining flatness once again. The analysis is based on the Malgrange preparation theorem. This approximation is valid in an entire neighborhood of a solution of the closure and shape matching constraints x 0 . The effect on ( 13) is the addition of a purely imaginary term

δ δx (S 1 (x) + S 2 (x) + S 3 (x) + iµ(x -x 0 ) 2 ) = θ x,1 + θ x,2 + θ x,3 + i 2µ(x -x 0 ) = 0 , ( 14
)
where µ depends on everything but x. To solve [START_REF] Dona | Numerical methods for eprl spin foam transitionamplitudes and lorentzian recoupling theory[END_REF] we need both the real and the imaginary part to vanish.

2 This is possible only if x is a solution of the constraints x = x 0 and the flatness condition is satisfied, θ x,1 + θ x,2 + θ x,3 = 0. Hellmann and Kamiński [START_REF] Hellmann | Holonomy spin foam models: Asymptotic geometry of the partition function[END_REF] suggested an argument for the flatness of the EPRL model based on a different strategy. First, they perform the summation over the bulk degrees of freedom and then they look for saddle points in the integrals over the holonomies. In the large boundary spins limit, the holonomies dominating the spin foam integrals are the ones contained in the wavefront set [START_REF] Hörmander | The Analysis of Linear Partial Differential Operators I: Distribution Theory and Fourier Analysis[END_REF] of the partition function. The wavefront set of W ∆3 fixes the product of the group elements around its bulk face g x = 1. If we parametrize g x = exp(iΘ m • σ/2) this condition imposes Θ = 0. Following [START_REF] Hellmann | Holonomy spin foam models: Asymptotic geometry of the partition function[END_REF] we can map3 g x into the geometrical holonomy around the triangle shared by the three 4-simplices. In this case, Θ can be interpreted as 2π minus the deficit angle around the internal face, obtaining the flatness condition.

Conclusion

We presented a numerical exploration of the semiclassical limit of the ∆ 3 transition amplitude in the 4D BF spin foam model. We computed the transition amplitude with coherent boundary data corresponding to Regge triangulations, both flat and curved. We found a stationary phase point in the sum over the bulk spin with value compatible with the area of the dual triangle of the prescribed triangulation in all the three cases.

This was at first surprising. The flatness problem argues that the partition function of the EPRL spin foam theory, in the semiclassical limit, is dominated by flat geometries. Usually, this is interpreted as an indication that the simplicity constraints, responsible for reducing a topological BF theory to General Relativity, are not imposed correctly. The same arguments, applied directly to SU(2) BF spin foam theory, indicate that curved geometries should be suppressed in the large boundary spins limit. This is in contradiction with our numerical result.

What is the resolution of this apparent tension? In Section 5 we analyzed the most common flatness arguments, and we discussed where they might be misleading. In particular, some of them derive their results from factoring many-vertices amplitudes as the product of single-vertex amplitudes and then applying the known asymptotics expansion [START_REF] Barrett | Asymptotic analysis of the eprl four-simplex amplitude[END_REF] to each vertex. They conclude that contributions of flat geometries dominate the stationary point analysis.

The same analysis can be applied to BF spin foam models as they have a similar asymptotic. The tension with our numerical calculation (which does not involve any approximation) must come from the saddle-point analysis of the many-vertices factorization. We suggest that the angles appearing in this analysis should not be interpreted as the deficit angles computed with a spin connection, i.e. as deficit angles à la Regge. This is indeed the simplest explanation for the case under examination. The equations of motion impose flatness of the BF connection but do not require the deficit angles computed with the piecewise flat metric (flat inside each 4-simplex) to vanish. The geometric deficit angles can assume any value, as indeed we observe in our numerical analysis.

In fact, the sums over the internal spins and the SU(2) integrals of the transition amplitude can be performed exactly, reducing it to a coherent {3nj} symbol. The large spin asymptotic of general coherent SU(2) invariants is well understood [START_REF] Dona | SU(2) graph invariants, regge actions and polytopes[END_REF][START_REF] Dona | Large spin limit of the EPRL-KKL spin foam vertex amplitude[END_REF]: the spin foam amplitude is not suppressed if the boundary data satisfy closure and shape matching constraints. The requirement of flat embedding in four-dimensional euclidean space is not necessary.

The same considerations may be partially applied to the EPRL spin foam model and variations, but the picture is less clear. Although these models have a semi-classical vertex expansion similar to BF models, the simplicity constraint are supposed to impose metricity on the spin connection. Care is needed in interpreting the angles appearing in the asymptotic analysis of EPRL amplitudes as geometric deficit angles. Our results suggest that this interpretation and the current flatness claims are misleading. To make any claim of flatness, we need to have a clear understanding of the correct semiclassical spin foam dynamics of many-vertices triangulations and its relation with the associated Regge equations of motions.

In a recent paper [START_REF] Asante | Effective Spin Foam Models for Four-Dimensional Quantum Gravity[END_REF], the authors explored the large spin limit of an effective spin foam model related to area-angle Regge calculus. They show that curved geometries are suppressed unless one considers a small Immirzi parameter in addition to large spins. They enrich their analysis with a numerical study confirming their findings. This calculation shows exactly the origin of the tension between our results and the various flatness arguments. These works share the common hypothesis that the asymptotic expression (11) holds also for extended triangulations. We have preliminary indications [START_REF] Dona | Asymptotic of spin foam amplitudes with many vertices[END_REF] that, while this expansion is correct for a single vertex, more care is needed in the presence of internal faces. All the flatness arguments become compatible with our numerical result if the angles in [START_REF] Magliaro | Emergence of gravity from spinfoams[END_REF] are not dihedral angles. What these arguments require to vanish is an angle related to the triangulation but that can not be interpreted as a pure deficit angle.

To conclude, our result is twofold. On the positive side, we showed that the flatness claims of the EPRL spin foam models are misleading. In the context of BF spin foam theory, we found an explicit counter-example in which they are not valid. We argue that the same should be valid for the EPRL model. We believe that flatness arguments do not invalidate the EPRL model as a good candidate for the LQG dynamics. On the negative side, we remark that the semiclassical limit of spin foam theories is not yet well understood. We find evidence that the geometric Regge action does not emerge from the large spin limit of SU(2) BF theory. Our result suggests that this could also be the case for other spin foam models like the EPRL model. Therefore we should not content ourselves with the present status of semiclassical spin foam analyses. We still lack a complete study showing that the true Regge action with the correct geometric interpretation arises for more involved triangulations. In our opinion, this is an urgent problem: a better understanding of the relation between quantum and classical discrete geometry will be one of the most important steps forward to make spin foam theories a valuable candidate for quantum gravity.

We believe that numerics will play a central role in future studies of spin foams. The numerical algorithm we described in this paper is powerful and can be applied to general amplitudes, with many vertices and internal faces, and different theories. We plan to use it to make a direct comparison of the large spin regime of BF and EPRL spin foam theories and finally provide an answer to this important question.

Addendum

After the publication of this paper, we continued developing our numerical code. The completely new version sl2cfoam-next [START_REF] Gozzini | High-performance spinfoam numerics[END_REF] improved on every aspect upon sl2cfoam, gaining several order of magnitudes in speed and efficiency. We expanded our case study, increasing the scale parameter to λ = 50 for the three boundary data studied in this paper. We included a third curved triangulation4 characterized by δ = 1.6 and scale λ = 50.

While the stationary phase point for flat boundary data is stable under rescaling, for curved boundary data we can infer a new trend. The characteristic oscillations of the partial sum around the "stationary phase point" can be many orders of magnitude greater than the value of the total amplitude, see Figure 6. This is evident for the new curved configuration we included here. We can interpret this as a signal of an exponential suppression of the full amplitude. We also studied the scaling behavior of the full amplitudes for the four boundary data, see Figure 7.

While the exponential suppression is particularly marked for the new curved boundary data, it is barely noticeable for the curved boundary data previously considered. In light of the recent results [START_REF] Engle | Comment on "EPRL/FK Asymptotics and the Flatness Problem[END_REF] we can interpret this difference as larger and smaller exponential suppression coefficients due to the accidental curvature constraint. Similar results can also be obtained for the EPRL model and will be published soon [START_REF] Gozzini | High-performance spinfoam numerics[END_REF].

We have to reconsider our conclusions. Our numerical analysis looks for the bulk spins that dominate the summation. We interpreted these points as stationary phase points. However, this interpretation seems to be inaccurate. We do not have numerical access to the action that enters the asymptotic analysis. We can infer the configurations that dominate the summation over bulk spins, but since these cannot be true stationary phase points in a discrete summation, we cannot derive from our analysis the exact scaling and the eventual exponential suppression of the resummed amplitude. We can find hints of the exponential suppression by noticing the presence of very large oscillations compared to the final value of the amplitude, but still the exponential suppression might be so slow that it could not be seen in the range of spins we have access to numerically. We conclude that, although still both flat and curved geometries dominate the sum over bulk degrees of freedom, the accidental curvature constraint in the curved case is present and results in an exponential suppression of the total amplitude in the large spin limit. From top left to bottom right the four plots correspond to boundary data with λ = 40 and δ = 0, δ = 3.60, δ = 2.47 and δ = 1.6. We also plot the power-law λ -12 for comparison. The scaling of the flat configuration is compatible with a power-law, while the scaling of the configuration with δ = 1.6 is exponentially suppressed. With the other two boundary data it is difficult to distinguish a power-law scaling from an exponential suppression.

Acknowledgments

The work of P.D. is partially supported by the grant 2018-190485 (5881) of the Foundational Questions Institute and the Fetzer Franklin Fund. We thank Lorenzo Bosi for his technical support with the CPT servers. F.G. acknowledges the Centre de Calcul Intensif d'Aix-Marseille for granting access to its high performance computing resources. We thank Carlo Rovelli for very useful discussions and for encouraging us in writing this paper. Special thanks go to Simone Speziale, for motivating us not to give up on this project, helping us with countless inputs and ideas, while taking care of another new little piece of his life. We thank Muxin Han, Hongguang Liu, and Jonathan Engle for suggesting us to expand our analysis and for pointing out the problem in our conclusions.

Appendices A Convention for the {15j} symbol

The {15j} symbol we use in this work is the irreducible {15j} symbol of first type, following the convention of [START_REF] Yutsis | Mathematical apparatus of the theory of angular momentum[END_REF]. The definition and its graphical representation is the following

j 1 j 5 k 1 j 3 k 4 l 5 l 1 l 4 l 2 l 3 k 1 j 2 k 5 j 4 k 3 - =    j 1 j 2 j 3 j 4 j 5 l 1 l 2 l 3 l 4 l 5 k 1 k 2 k 3 k 4 k 5    = (15)
= (-1)

5 i=1 ji+li+ki s (2s + 1) j 1 k 1 s k 2 j 2 l 1 j 2 k 2 s k 3 j 3 l 2 (16) × j 3 k 3 s k 4 j 4 l 3 j 4 k 4 s k 5 j 5 l 4 j 5 k 5 s j 1 k 1 l 5 . (17)

B Coherent states

Coherent states are a fundamental ingredient to study the semiclassical behavior of a spin foam amplitude. We briefly review the conventions used in this paper. For more details we refer to [START_REF] Perelomov | Generalized coherent states and their applications[END_REF] or to other numerical works [START_REF] Dona | SU(2) graph invariants, regge actions and polytopes[END_REF] and [START_REF] Dona | Numerical methods for eprl spin foam transitionamplitudes and lorentzian recoupling theory[END_REF]. A SU(2) coherent state |j, n in the irreducible representation of spin j is given by the action on the lowest weight state |j, -j of a SU(2) group element g n corresponding to the rotation that transforms the unitary vector z = (1, 0, 0) into n = (sin Θ cos Φ, sin Θ sin Φ, cos Θ). This coherent state minimizes the uncertainty in the two directions orthogonal to n. The Euler angles of group element g n are related to the angles Θ and Φ by

g n = e -iΦ σz 2 e -iΘ σy 2 e iΦ σz 2 , (18)
where we expressed the SU(2) generators i σ/2 in terms of the Pauli matrices. The scalar product between a coherent state and a state in the magnetic basis is given by j, m| j, n = j, m| D

(j) (g n ) |j, n = D (j) m-j (g n ) = D (j) m-j (Φ, Θ, -Φ) , (19)
where D

Finally we can define a coherent tetrahedron as the projection on the intertwiner space of the tensor product of four SU(2) coherent states. Each one of them is associated to a face of the tetrahedron with area j i and normal n i c j12 ( n i ) : = j i ; j 12 | j 1 , n 1 ; j 2 , n 2 ; j 3 , n 3 ; j 4 , n 4 (

) 24
The graphical representation of an SU(2) coherent state is given by a line ending on a little circle while the coherent states coefficients are represented by (

|j, n = c j12 ( n i ) = j 1
) 28
To keep the picture as simple as possible we opted for not using the graphical representation of the coherent tetrahedra c ia ( n f a ) that we factorize. The "dipole" like diagrams contribute with a phase and the inverse dimension of the intertwiner space [START_REF] Gozzini | High-performance spinfoam numerics[END_REF] with χ = 2(j 123 + j 234 + j 124 + j 134 + j 456 + j 156 + j 346 + j 356 + j 235 ) + j 123 + j 345 + j 156 . This phase is obtained from the contributions of ( 29) and (-1) 2j required to invert the orientation of an internal line of spin j. We finally derive [START_REF] Conrady | On the semiclassical limit of 4D spin foam models[END_REF].

Figure 1 :

 1 Figure 1: Left: The 2-complex of the ∆3 triangulation. We highlight in red the bulk face. Right: The ∆3 spin foam amplitude in graphical notation. We are omitting a summation over the internal spin x and a dimensional factor 2x + 1. The boxes represent integrals over the SU (2) group, and the strands are irreducible representations labeled by a spin j abc . The empty dots on the boundary represents the SU(2) coherent states.

  i a are the boundary intertwiners, k b are the internal intertwiners and the sum over the spin x is bounded by triangular inequalities. The complex coefficients c ia ( n f a ) are the Livine-Speziale coherent intertwiners in the recoupling basis, they depend on the normals n f a and areas of the tetrahedron a (see appendix B for their definition).

Figure 2 :

 2 Figure 2: Left: The geometry of the ∆3 triangulation. The numbered circles correspond to points. Lines correspond to segments. Each color corresponds to a different 4-simplex. The bulk triangle is highlighted in red. Right: The three 4-simplices are shown separately.

Figure 5 :

 5 Figure5: Second curved configuration. Left: Averaged running sum Rw(x) with interval sizes 2,6, and 10. We normalize respect to the value of its maximum. Right: Real part of the partial sum Pw(x) normalized respect the absolute value of the amplitude. In both plots we mark with a red line the value of the geometrical area of the bulk triangle, and with a green band the values of the numerical estimate of the stationary phase point within one standard deviation.

Figure 6 :

 6 Figure 6: Real part of the partial sum Pw(x) normalized with respect to the absolute value of the amplitude. In all plots we mark with a red line the value of the geometrical area of the bulk triangle. From top left to bottom right the four plots correspond to boundary data with λ = 50 and δ = 0, δ = 3.60, δ = 2.47 and δ = 1.6.

Figure 7 :

 7 Figure7: Log-Log plot of the absolute value of the W∆ 3 amplitude as a function of the scale λ. From top left to bottom right the four plots correspond to boundary data with λ = 40 and δ = 0, δ = 3.60, δ = 2.47 and δ = 1.6. We also plot the power-law λ -12 for comparison. The scaling of the flat configuration is compatible with a power-law, while the scaling of the configuration with δ = 1.6 is exponentially suppressed. With the other two boundary data it is difficult to distinguish a power-law scaling from an exponential suppression.

2 dj 1 j 2 j 3 j 4 m 1 m 2 m 3 m 4 (j 1 j 2 j 12 m 1 m 2 m j 12 j 3 j 4 -m m 3 m 4 .

 2444 (Φ, Θ, -Φ) is the Wigner matrix element of g n parametrized in terms of its Euler angles. These conventions are used by sl2cfoam and Wolfram Mathematica. SU(2) coherent states form an overcomplete basis in the irreducible representation of spin j space1 (j) = (2j + 1) S n |j, n j, n| ,(20)where d n is the normalized measure on the sphere. A 4-valent intertwiner in the recoupling basis is given by|j i ; j 12 = m1 j12) |j 1 m 1 |j 2 m 2 |j 3 m 3 |j 4 m 4(21)and the (4jm) symbol is given by the contraction

  m1-j1 (g n1 )D (j2) m2-j2 (g n2 )D (j3) m3-j3 (g n3 )D (j4) m4-j4 (g n4 ) .

C(

  Derivation of the ∆ 3 spin foam amplitudeIn this appendix we derive the formula (5) for the spin foam transition amplitude associated to the ∆ triangulation starting from its integral representation. We use the SU(2) graphical calculus to represent the amplitude asW ∆3 (j f , n f a ) =x we picked a conventional orientation of the faces that are capped by coherent states on both ends, represented by empty circles. We label the 4-simplices using the five points, the edges, dual to a tetrahedron, are labeled by four points and the faces, dual to triangles, are labeled by three points. We highlight in red the internal face x.We perform the integrations over SU (2) and we get W ∆3 (j f , n f a ) =x 2i a + 1)c ia ( n f a ) ×

  invert the orientation of some internal lines to conform with[START_REF] Speziale | Boosting Wigner's nj-symbols[END_REF] and obtain the ∆ 3 transition amplitude in terms of three {15j} symbolsW ∆3 (j f , n f a ) = (-1) χ x (-1) x (2x + 1) k b b (-1) k b (2k b + 1)ia a (2i a + 1)c ia ( n f a )

345 k 1 j 145 i 2 j 235 x j 134 j 124 j 245 j 125 k 2 j 123 i 3 j 234

	 	 	i 4 j
		

123 k 2 j 235 i 5 j 136 x j 125 j 256 j 236 j 356 k 3 j 156 i 6 j 126

  123 +j 234 +j 124 +j 134 +j 456 +j 156 +j 346 +j 356 +j 235 )+j 123 +j 345 +j 156 is a phase function only of the boundary spins,

				
				
				
	×	 	i 7 j 156 k 3 j 136 i 8 j 145 x j 356 j 346 j 146	 
			j 134 k 1 j 345 i 9 j 456	
	where χ = 2(j

Table 1 :

 1 We compare the numerical estimate of the stationary phase point with the geometrical area at different scales.

		1	13.4	2	0.2
	λ = 12 :	16 ± 1	16.1	6	0.4
	λ = 14 :	18 ± 1	18.8	16	0.9
	λ = 16 :	21 ± 1	21.5	33	3
	λ = 18 :	23 ± 1	24.1	65	7.8
	λ = 20 :	26 ± 1	26.8	121	21
	λ = 30 :	40 ± 1	40.2	∼ 2500	∼ 2700

The ∆ 3 triangulation is formed by three 4-simplices, sharing a common triangle. The boundary of the triangulation consists of nine tetrahedra and the bulk consists of three tetrahedra all sharing the common triangle x.

The same conclusion is obtained in[START_REF] Engle | private communications and ILQGS seminar[END_REF] in a more rigorous way. Instead of looking for stationary phase points they perform the summation over x using the Poisson resummation. The flatness condition is derived by requiring a non exponentially suppressed sum.

For Regge boundary data, the coherent {15j} symbol has two distinct saddle points. Up to gauge, the product of the two group elements on a wedge at the saddle point is characterized by ± the dihedral angle associated with that wedge. Each solution is responsible for one exponential forming the cosine in[START_REF] Magliaro | Emergence of gravity from spinfoams[END_REF]. If we take the saddle points of the three vertices with angles with the same sign, the angle of gx is ± the sum of dihedral angles around the shared triangle. With these two group elements, we use the map between SU (2) × SU (2) and SO(4) to derive a non trivial SO(4) geometrical holonomy, that is a rotation of angle Θ on the plane orthogonal to the shared triangle.

this case was suggested to us by Hongguang Liu in a discussion.