Pietro Donà

Marco Fanizza

Giorgio Sarno

Simone Speziale

Numerical study of the Lorentzian Engle-Pereira-Rovelli-Livine spin foam amplitude

The Lorentzian Engle-Pereira-Rovelli-Livine spin foam amplitude for loop quantum gravity is a multidimensional noncompact integral of highly oscillating functions. Using a method based on the decomposition of Clebsch-Gordan coefficients for the unitary infinite-dimensional representations of SL(2,C) in terms of those of SU(2), we are able to provide for the first time numerical evaluations of the vertex amplitude. The values obtained support the asymptotic formula obtained by Barrett and collaborators with a saddle point approximation, showing, in particular, a power law decay and oscillations related to the Regge action. The comparison offers a test of the efficiency of the method. Truncating the decomposition to the first few terms provides a qualitative matching of the power law decay and oscillations. For vector and Euclidean Regge boundary data, a qualitative matching is obtained with just the first term, which corresponds to the simplified EPRL model. We comment on future developments for the numerics and extension to higher vertices. We complete our work with some analytic results: We provide an algorithm and explicit configurations for the different geometries that can arise as boundary data, and explain the geometric consequences of the decomposition used.

Introduction

We report on numerical studies of the Lorentzian Engle-Pereira-Rovelli-Livine (EPRL) vertex amplitude for spin foam quantum gravity [START_REF] Engle | LQG vertex with finite Immirzi parameter[END_REF][START_REF] Freidel | A New Spin Foam Model for 4D Gravity[END_REF] (for reviews, see [START_REF] Perez | The Spin Foam Approach to Quantum Gravity[END_REF][START_REF] Rovelli | Covariant Loop Quantum Gravity[END_REF]). These are aimed in particular at testing its large spin asymptotic behavior, as derived by Barrett et al. [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF]. The asymptotics can be interpreted as a semiclassical limit for the vertex amplitude and shows the presence of the Regge action, for a suitable choice of boundary data. It establishes an important relation to (discrete) general relativity that is at the heart of most of the physical applications of the EPRL model, thereby the interest in numerical verifications.

Our numerical work is based on the code sl2cfoam developed in [START_REF] Dona | Numerical methods for EPRL spin foam transition amplitudes and Lorentzian recoupling theory[END_REF] and publicly available at [START_REF]The code sl2cfoam is public and available at the[END_REF], optimized here for the specifics of the vertex amplitude. The code implements the method of [START_REF] Speziale | Boosting Wigner's nj-symbols[END_REF] to compute Lorentzian Clebsch-Gordan coefficients and spin foam amplitudes. Our goal is twofold: on the one hand, to provide the first numerical test of Barrett's Lorentzian asymptotics results [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF], thus verifying the saddle point approximation underpinning it, as well as estimating its validity. On the other hand, to examine the robustness of sl2cfoam and test its performances. This means being able to compute Lorentzian spin foam amplitudes in the deep quantum regime, a necessary step for the spin foam formalism.

The results we found are encouraging, albeit not completely satisfactory. The numerics confirm Barrett's formula for boundary data corresponding to vector geometries or to Euclidean 4-simplices. In these cases, the asymptotics is reached fast, with an error below 10% at spins of order 10. This confirms the validity of the saddle point approximation already at small spins, a situation similar to the much more extensively studied case of asymptotics of SU(2) invariants, see e.g. [START_REF] Ponzano | Semiclassical limit of Racah coefficients[END_REF][START_REF] Donà | SU(2) graph invariants, Regge actions and polytopes[END_REF]. In parallel, it shows the accuracy and robustness of sl2cfoam. To give an idea of the heavy computational load required to achieve these results, the plot used to confirm the Euclidean asymptotics at spins of order 10 needed the evaluation of 10 11 configurations, each of which is determined by up to 10 4 one-dimensional unbounded integrals of oscillating functions. We used a server provided by the CPT with 32 Intel R Xeon R CPUs E5-2687W v2 at 3.40GHz with 264 gigabytes of RAM. The evaluation times depend on the boundary data, as we will explain in the next sections. For an Euclidean 4-simplex, the main plot of section 5.3 took approximately 3 weeks. For vector geometries the two plots in section 5.2 took about 30h each. Lorentzian 4-simplices are even more demanding, and a bottleneck in our computational precision restricted the number of data points accessible to us. These are shown in section 5.4, which took about 2 days. These were enough to provide evidence of a critical behavior with the predicted power law falloff, as well as nonmonotonic behavior depending on the Immirzi parameter. But a more significative agreement with the asymptotic formula, in particular seeing the predicted oscillations, was not obtained, and will require more computational work. Our results also show that the simplified model introduced in [START_REF] Speziale | Boosting Wigner's nj-symbols[END_REF] captures the key features of the full EPRL asymptotics for vector and Euclidean Regge geometries, but fails to do so, and instead decays exponentially, for Lorentzian geometries.

The increased difficulty of Lorentzian 4-simplices has partially to do with the fact that one cannot resort to an equispin configuration. The best we found was a certain isosceles configuration, but even then one has to deal with individual spins that are much higher than the value of the rescaling parameter controlling the asymptotics. Furthermore, the internal sums introduced by the method of [START_REF] Speziale | Boosting Wigner's nj-symbols[END_REF] turn out to converge more slowly than for a Euclidean 4-simplex. This fact actually has a geometric origin, which explains how the Lorentzian 4-simplex is reconstructed from many Euclidean 4-simplices compatible with the SU(2) {15j} symbol, and whose tetrahedra are suitably transformed thanks to the half-edge booster functions. This discussion also clarifies the critical behavior observed for the simplified model.

To obtain the numerical results presented, it proved necessary to develop an efficient algorithm to construct the boundary data. This is a spin-off of our work, and we devote a section of the paper to explain how this problem can be solved and provide in the same online repository [START_REF]The code sl2cfoam is public and available at the[END_REF] the relevant pieces of code. For the vector geometries, we used the explicit parametrization worked out in [START_REF] Donà | SU(2) graph invariants, Regge actions and polytopes[END_REF]. For the Regge geometries, we used a "deconstruction" algorithm, whereby we start from a geometric 4-simplex identified by its vertices, and derive the 3D data associated with it by the EPRL model's Y map. This procedure exposes some useful properties of the 3D data and was helpful for us to better understand the details of the saddle point analysis of [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF]. We provide some more explicit formulas for the Hessian at the saddle point, needed for the numerical comparison of the analytic formula with the data points.

Barring the increased computational power to confirm precisely the analytic formula for Lorentzian 4simplices, our work confirms the validity of this crucial result for the semiclassical limit of the EPRL model, but more importantly supports the robustness of the method [START_REF] Speziale | Boosting Wigner's nj-symbols[END_REF] and code sl2cfoam [START_REF] Dona | Numerical methods for EPRL spin foam transition amplitudes and Lorentzian recoupling theory[END_REF], which can now be used to perform more ambitious calculations. To that end, we point out the heavy cost of using coherent intertwiners. These are needed for the asymptotic formula, but many physical applications of the EPRL spin foam model can be done in the much cheaper orthonormal basis of intertwiners. Some of us presented calculations with two nonsimplicial vertices in [START_REF] Sarno | 2-vertex Lorentzian Spin Foam Amplitudes for Dipole Transitions[END_REF]. The present code is much more performing than the one there used, and we hope to apply it in future work to calculations with a few simplicial vertices and orthonormal intertwiners.

With the intent to make the paper concise and its technical results easier to appreciate, we refer to the existing literature for the necessary background. In particular, to [START_REF] Engle | LQG vertex with finite Immirzi parameter[END_REF][START_REF] Livine | Consistently Solving the Simplicity Constraints for Spinfoam Quantum Gravity[END_REF][START_REF] Freidel | A New Spin Foam Model for 4D Gravity[END_REF] for the analytic background on the EPRL model, to [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF] for its saddle point approximation (see also [START_REF] Barrett | A Lorentzian signature model for quantum general relativity[END_REF][START_REF] Barrett | Asymptotic analysis of the EPRL four-simplex amplitude[END_REF][START_REF] Conrady | On the semiclassical limit of 4D spin foam models[END_REF][START_REF] Barrett | Quantum gravity asymptotics from the SU(2) 15j symbol[END_REF][START_REF] Han | Asymptotics of Spinfoam Amplitude on Simplicial Manifold: Lorentzian Theory[END_REF] for related work), to [START_REF] Speziale | Boosting Wigner's nj-symbols[END_REF] for the method used for the exact evaluation, and to [START_REF] Dona | Numerical methods for EPRL spin foam transition amplitudes and Lorentzian recoupling theory[END_REF] for the numerical code, describing its details, functioning scheme, and generic performance tests. Section 2 contains the minimal information about the vertex amplitude for the paper to be self-contained and gives an overview of the approach used for the calculation. Section 3 describes our algorithm to identify boundary data with the required properties, and the explicit configurations used in the numerical calculations. Section 4 discusses the central point of the convergences of the internal sums and presents a useful approximation that we use to shorten significantly the numerical time of some numerical evaluations. Section 5 contains our main results, a selection of data and their comparison with the analytic asymptotic behavior. Section 6 contains a discussion on the role of the booster functions in mapping Lorentzian 4-simplices from the Euclidean ones associated with the SU(2) {15j} symbol, as well as a summary of the numerical situation and future developments. After our conclusions, we supplement the paper with two Appendices. In one we provide explicit formulas to fix our conventions. In the second we summarize the asymptotic analysis of the EPRL vertex [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF] and report on the form of the Hessian at the critical point.

EPRL vertex amplitude and its asymptotic limit

The (coherent) vertex amplitude for the Lorentzian EPRL spin foam model [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 SL(2, C) invariant associated to the boundary graph of a 4-simplex, A v (j ab , n ab ) := 5 a=2 dh a a<b D (γj ab ,j ab ) j ab ,-n ab ,j ab , n ba (h -1 a h b ), [START_REF] Engle | LQG vertex with finite Immirzi parameter[END_REF] where h a ∈ SL(2, C), and D (ρ,k) (h) are infinite-dimensional unitary representations of the principal series, labeled by ρ ∈ R and k ∈ Z/2. The irreducible representations (irreps) are restricted to satisfy ρ ab = γk ab , where γ is the Immirzi parameter, and expressed using Naimark's canonical basis with minimal spin eigenvalues j ab = k ab , and unit vectors n ab used to build coherent intertwiners [START_REF] Livine | A New spinfoam vertex for quantum gravity[END_REF][START_REF] Livine | Consistently Solving the Simplicity Constraints for Spinfoam Quantum Gravity[END_REF]. Conventions and explicit formulas are reported in Appendix A. The minus sign in the vectors appearing as rows is for later convenience. The indices a, b = 1, . . . 5 stand for the nodes of the graph, and notice that one redundant integration has been removed, which is necessary to ensure finiteness [START_REF] Engle | Regularization and finiteness of the Lorentzian LQG vertices[END_REF][START_REF] Baez | Integrability for relativistic spin networks[END_REF] (see also the discussion in [START_REF] Sarno | 2-vertex Lorentzian Spin Foam Amplitudes for Dipole Transitions[END_REF]).

The boundary data of the amplitude are ten spins j ab = j ba and twenty unit vectors n ab = n ba . Four different subsets of data play an important role, see Fig. 1: First, closed twisted geometries, defined by the following closure conditions, b =a j ab n ab = 0, ∀a = 1.

(

) 2
This means that the spins and normals data around each node define a tetrahedron. 1Within the space of closed twisted geometries, we can distinguish two disconnected subsets. One is vector geometries, if furthermore there exist five rotations R a ∈ SO(3) such that the two normals associated to each face can be made opposite to one another:

R a n ab = -R b n ba . (3)
These rotations are determined by h a at the corresponding critical point. Clearly, the simplest possible case is when the normals are pairwise antiparallel,

n ab = -n ba , (4)
which plays a useful role in choosing data. It is a remarkable fact that all Euclidean Regge geometries, namely all Euclidean 4-simplices in this context, are a subset of vector geometries. They are characterized by additional shape-matching constraints on top of Eq. ( 3), which can be described either using bivectors [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF][START_REF] Barrett | Relativistic spin networks and quantum gravity[END_REF] or spherical cosine laws [START_REF] Dittrich | Area-angle variables for general relativity[END_REF][START_REF] Donà | SU(2) graph invariants, Regge actions and polytopes[END_REF]. Their neat implication is to guarantee that the five tetrahedra determined by the closure conditions coincide with those determined as the boundary of a 4-simplex having the spins as areas. 2 What happens in the simplest case (4) is that there is always one critical point with all holonomies at the identity, and only when the shape-matching conditions are further satisfied does one find a second non trivial critical point that allows the reconstruction of the Euclidean 4-simplex.

The second disconnected subset of the closed twisted geometries are the Lorentzian Regge geometries, for which the boundary data coincide with those determined by a Lorentzian 4-simplex with all tetrahedra spacelike. The subset can again be characterized by shape-matching constraints in terms of bivectors [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF][START_REF] Barrett | Relativistic spin networks and quantum gravity[END_REF] or Lorentzian spherical cosine laws [START_REF] Anzà | A note on the secondary simplicity constraints in loop quantum gravity[END_REF][START_REF] Dona | Asymptotics of lowest unitary SL(2,C) invariants on graphs[END_REF]. Furthermore, while in this subset [START_REF] Perez | The Spin Foam Approach to Quantum Gravity[END_REF] are not satisfied, it is still true that the 3D normals are related by five complex transformation H a such that

H a n ab = -H b n ba . (5)
These unimodular complex matrices are non unitary 3D representations of Lorentz boosts. They are again determined by the h a at the corresponding critical point. How these conditions allow the embedding of the boundary data in a flat Euclidean or Lorentzian 4simplex will be explained in detail below. It is in the meantime useful to remark that the areas alone suffice to characterize a 4-simplex and distinguish whether it is Euclidean or Lorentzian. Data with areas compatible with a Euclidean 4-simplex can be either a vector geometry or a Euclidean Regge geometry depending on the compatibility or lack thereof of the normals with the Euclidean 4-simplex. This requirement can be expressed in terms of the shape-matching conditions reducing a twisted geometry to a Regge geometry [START_REF] Dittrich | Area-angle variables for general relativity[END_REF][START_REF] Donà | SU(2) graph invariants, Regge actions and polytopes[END_REF].

Generic data namely open twisted geometries 35d

The coherent amplitude ( 1) is a complex number, and it is important to make a remark about its phase before continuing. The original formulation of the model [START_REF] Engle | LQG vertex with finite Immirzi parameter[END_REF] uses the standard basis of orthonormal intertwiners, which are labeled by a spin number i a per vertex, instead of the coherent intertwiners, labeled by unit vectors n ab . The resulting vertex amplitude A v (j ab , i a ) is a real number, provided a certain phase convention for the Lorentz irrep matrices is chosen [START_REF] Speziale | Boosting Wigner's nj-symbols[END_REF]. The interest in using the complex coherent amplitude (1) comes from the finer characterization of the boundary data and the clearer geometric interpretation that it permits. However, while the original amplitude is SU(2) invariant, Eq. ( 1) is only SU(2) covariant: A rotation of the 4-normals on a node preserves the amplitude only up to a phase. 3 This phase is irrelevant for geometric considerations because the coherent states themselves are only defined up to a phase, thus for any chosen boundary configuration, one can suitably tune the coherent state phase to have a real amplitude. Accordingly, we can neglect this phase and factor out the rotational gauge freedom at each node in our geometric considerations. Up to rotations, the boundary data span a 35-dimensional space. The closed twisted geometries are a 20-dimensional subset, 10 areas and 2 × 5 = 10 shape parameters, the vector geometries a 15-dimensional subset, ten areas and five angles [START_REF] Donà | SU(2) graph invariants, Regge actions and polytopes[END_REF], and the Regge geometries a ten dimensional subset, ten lengths or ten areas (up to some discrete ambiguity in the Euclidean case).

Thanks to properties of the coherent states, integrands like Eq. ( 1) can be written in exponential form with a linear dependence on the spins [START_REF] Livine | A New spinfoam vertex for quantum gravity[END_REF]. In particular, under a homogeneous rescaling j ab → λj ab we can write the integrand as exp λS(g a ; j ab , n ab ). This gives the possibility of evaluating the amplitude with a saddle point approximation in the homogeneous large spin limit, λ → ∞. This approximation was studied in [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF], 4 and it was found that the leading order asymptotics of the coherent amplitude (1) has the following dependence on the boundary data:

1. For generic open or closed twisted geometries: exponential decay in λ; 2. For vector geometries: power law decay λ -12 ; 3. For Regge geometries: power law decay λ -12 with nonambiguous real oscillations, and frequency given by the Regge action. The Regge action is independent of γ for Euclidean data and depends linearly on it for Lorentzian data.

Our goal in this paper is to test numerically these saddle point approximations, in particular the power law decays and the frequency of the oscillations. Such a numerical evaluation is not an easy task: The formula for the amplitude contains 12 unbounded integrations, and a complex and rapidly oscillating integrand built out of sums of products of hypergeometric functions. A direct approach using for instance adaptive Monte Carlo methods5 appears daunting, and we are not aware of any results even for the much simpler Lorentzian Barrett-Crane integrals. Facing these difficulties, the idea of [START_REF] Speziale | Boosting Wigner's nj-symbols[END_REF] was to reduce as much as possible the unbounded integrations, taking advantage of the known factorization of SL(2, C) Clebsch-Gordan coefficients into SU(2) ones. Applied to Eq. ( 1), the method in [START_REF] Speziale | Boosting Wigner's nj-symbols[END_REF] splits the expression into a convolution between an SU(2) {15j} symbol and one-dimensional boost integrals, called booster functions and denoted B γ 4 :

A v (j ab , n ab ) = l ab ≥j ab ka,ia

Here d j := 2j + 1. The {15j} symbol is labeled by internal spins l ab and intertwiners k a , except at the node 1 without group integration, where it is labeled directly by the boundary spins and intertwiners. The unbounded sums over l ab go from j ab to infinity, whereas the sums over the intertwiner labels k a , i a are bounded by the usual SU(2) triangle inequalities. We define the simplified EPRL model consistently with [START_REF] Speziale | Boosting Wigner's nj-symbols[END_REF][START_REF] Dona | Numerical methods for EPRL spin foam transition amplitudes and Lorentzian recoupling theory[END_REF][START_REF] Sarno | 2-vertex Lorentzian Spin Foam Amplitudes for Dipole Transitions[END_REF] as the truncation of the EPRL amplitude [START_REF]The code sl2cfoam is public and available at the[END_REF] where only the first term in the summation l ab = j ab is considered.

The booster functions depend both on the internal and the boundary data, and to obtain the coherent amplitude (1), they are contracted with the coefficients c i of the coherent intertwiners determined by the normals n ab . The graphical notation used in the second line helps to keep track of which links carry internal spins, as well as of the orientation of the normals. The booster functions encode the details of the Y -map defining the EPRL model, and are defined as follows,

B γ 4 (j f , l f ; i, k) = d i d k m f j f m f (i) l f m f (k) 1 4π ∞ 0 dr sinh 2 r 4 f =1 d (γj f ,j f ) j f l f m f (r) = j 4 i j 3 j 2 j 1 l 4 k l 3 l 2 l 1 , (8)
in terms of generalized four-legged (4jm) Wigner symbols and boost matrix elements d (ρ,k) (r), see [START_REF] Speziale | Boosting Wigner's nj-symbols[END_REF] for explicit formulas. To not be confused with d i and d k the dimension of the SU(2) irrep of spin i and k, see Appendix A for more details on our conventions. For all links outgoing, the coherent intertwiner coefficients are

c i ( n f ) := m f j f m f (i) 4 f =1 D (j f ) m f j f ( n f ) = ñ1 ñ2 ñ3 ñ4 j 1 j 2 j 3 j 4 i , (9)
where

D (j) ( n) = e -iΦJz e -iΘJy e iΦJz (10)
is a Wigner matrix with Θ and Φ the zenithal and azimuthal angles of n. If one of the links is incoming one multiplies the above formula by (j) mn = (-1) j-m δ m,-n . In spite of its complexity, the vertex amplitude ( 6) is in a form suitable for numerical evaluations. The unbounded integrations have been limited to the booster functions, that are accessible numerically, but at the price of introducing infinite summations, the sums over the internal spins l, or "internal sums" for short. These may be just as hard to deal with numerically if they cannot be efficiently truncated. We will discuss this point below in detail. Let us first examine the different ingredients required in setting up the algorithm.

The first is evaluating the {15j}. Optimized algorithms exist in the literature for the simpler {6j} and {9j} symbols, but the {15j} can be significantly more complicated. In fact, the problem would be intractable if one had to compute it starting from its definition as ten (3jm) symbols, contracted following the graph drawn in Eg. [START_REF]The code sl2cfoam is public and available at the[END_REF]. We can instead exploit the invariance of the sum over intertwiners under recoupling moves to trade the irreducible {15j} appearing in Eq. ( 7) with a reducible one, a strategy successfully applied in [START_REF] Donà | SU(2) graph invariants, Regge actions and polytopes[END_REF] to compute numerically the semiclassical limit of the coherent {15j} symbol. The simplest choice is a reducible {15j} given by one sum of five {6j} symbols. The {6j} symbols can then be very efficiently computed thanks to the libraries wigxjpf developed in [START_REF] Johansson | Fast and accurate evaluation of wigner 3j, 6j, and 9j symbols using prime factorisation and multi-word integer arithmetic[END_REF], based on smart storage of factorials used to compute the symbols. With this procedure, computing a single {15j} for spins of order 10 takes about 10 -6 s. We refer to the original paper for a detailed discussion of the library's performances, accuracy, and memory management.

The second is the booster functions. Here and in sl2cfoam we use the definition (8) with a polynomial expression for the boost matrix elements in e -r derived in [START_REF] Collet | A (simple) expression of the unitary-irreducible SL(2,C) representations as a finite sum of exponentials[END_REF], which is more efficient to implement numerically than the usual expression as hypergeometric functions. 6 This expression is then amended with the phase conventions of [START_REF] Speziale | Boosting Wigner's nj-symbols[END_REF] to make the booster functions and vertex amplitudes with orthonormal intertwiners real. We compute the numerical integrals at fixed magnetic labels m f , then resum over the (4jm) coefficients. This gives an upper estimate of 4λ 3 integrals, since only three magnetic labels are linearly independent. The results are then stored using a hash table. With this procedure, computing a booster for spins of order 10 takes approximately 6 minutes on a single core, and needs to be computed only once. The main limitation of this method is that numerical instabilities become important at spins of order 50. These instabilities prevent the convergence of the boost integrals and are presumably due to the algebraic manipulations of very large numbers in the integrand. We worked with the GNU MPFR library for multiple-precision floating-point computations. See [START_REF] Dona | Numerical methods for EPRL spin foam transition amplitudes and Lorentzian recoupling theory[END_REF] for explicit formulas and more details. This makes the boosters the most delicate part of our code, and an obvious target for improvements, as we will discuss below in Sec. 6.2.

The third is the coherent intertwiner coefficients [START_REF] Ponzano | Semiclassical limit of Racah coefficients[END_REF]. This is the simplest part of the code, obtained contracting the Wigner matrices [START_REF] Donà | SU(2) graph invariants, Regge actions and polytopes[END_REF] with (4jm) symbols computed with wigxjpf. The evaluation time of a single coherent state coefficient in sl2cfoam , is negligible with respect to the other computation times.

Finally one has to tabulate all these quantities and combine them to compute the amplitude. This is a delicate task that is performed using the hash-table algorithm described in detail in [START_REF] Dona | Numerical methods for EPRL spin foam transition amplitudes and Lorentzian recoupling theory[END_REF]. The resulting code is available in [START_REF]The code sl2cfoam is public and available at the[END_REF], and can be downloaded and used to independently rederive the results presented here, or push them forward using more powerful computers.

Constructing the boundary data

Our code allows us to compute the vertex amplitude for any boundary data, in principle. However, some choices are more efficient than others: for instance, symmetric configurations maximize the number of data points with a given maximal value of the rescaling parameter λ, thus optimizing the numerical effort. This is simply a consequence of the discreteness of the areas, and the fact that the numerical costs are strongly determined by the highest spin in the amplitude. For twisted, vector and Euclidean Regge geometries, one has access to the most symmetric configuration with all equal spins j ab ≡ j, and this is the one we focused on.

For Lorentzian Regge geometries, it is not possible to take all spins as equal since an equilateral Lorentzian 4-simplex does not exist. We then devised a Mathematica code to span Lorentzian 4-simplices with all tetrahedra spacelike, and we scanned the space of configurations with integer or half-integer areas. We found that the case with minimal highest spin is a configuration with six spins equal to 2, and four spins equal to 5, which can then be rescaled by an integer multiple (but not a half-integer). The four highest spins can all be associated to the same tetrahedron, and it is then best to choose to gauge fix the SL(2, C) integration on the same tetrahedron, so to avoid having to compute the heaviest booster function.

In all cases, we define the homogeneous rescaling parameter λ as an integer multiplying the smallest allowed configuration, hence j ab = λ/2 for equispin data, and j ab = (2λ, 5λ) for Lorentzian boundary data. The two spin configurations are summarized in Table 1. From a purely numerical perspective, these configurations are rather generic, in the sense that they test generic structures of the numerical code.

Twisted, vector, Euclidean data All spins:

λ/2 = 1 2 , 1, 3 2 , 2, . . .

Lorentzian data

Four spins in gauge-fixed tetrahedron: 5λ = 5, 10, . . .

Remaining six spins: 2λ = 2, 4, . . .

Table 1:

The two configurations used for boundary Regge data throughout the paper. The equispin data used for the Euclidean Regge geometries can also be used for twisted and vector geometry data: It is the choice of normals that distinguish the three cases.

The information just given suffices to follow the numerical analysis in the paper. On the other hand, the actual construction of the boundary data turned out to be a nontrivial part of our analysis, and for the interested reader, we explain in the rest of this section how it can be done, and provide more details on the specific configurations we chose. The online repository [START_REF]The code sl2cfoam is public and available at the[END_REF] contains the Mathematica notebooks that we used to construct the boundary data we used in our asymptotic analysis. The notebook can be used to generate the boundary data for any spin foam vertex amplitude.

The first thing to keep in mind is that the amplitude's boundary data (j ab , n ab ) are purely 3-dimensional: j ab label irreps of the canonical SU(2) matrix subgroup corresponding to rotations preserving the time direction t I := (1, 0, 0, 0), and n ab are unit vectors in R 3 . We then have both an algebraic and a geometric mapping of the data into four-dimensional structures. The algebraic mapping is straightforward: The stabilizer of t I is also the canonical SU(2) subgroup of SL(2, C) used in Naimark's basis of the unitary representation of the principal series (ρ, k) used in the EPRL model. Then the mapping is given by the map Y : j → (γj, j), and |j, n → |γj, j; j, n . See [START_REF] Engle | LQG vertex with finite Immirzi parameter[END_REF][START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF][START_REF] Martin-Dussaud | A Primer of Group Theory for Loop Quantum Gravity and Spin-foams[END_REF] for more details. The geometric mapping is based on looking at the 3D geometry described by the data, and studying when and how it can be taken as the boundary of a flat 4-simplex. This mapping is finer and allows the classification of the data reported in Fig. 1. Given a set of data (j ab , n ab ), it can be directly tested whether it describes one of the special geometric classes. For our purposes though, we are interested in the reverse, "deconstructing" problem: start from a geometric object of interest, like a Euclidean or Lorentzian 4-simplex, and deduce the 3D data to be used as boundary states.

A closed twisted geometry

To generate the data of a closed twisted geometry it is convenient to start with five tetrahedra with areas matching along the graph of the 4-simplex. Each one can be characterized by the four areas and two independent 3D dihedral angles, φ e ≡ φ (a) bc at the edge e of the tetrahedron a identified as the one shared with the tetrahedra b and c in the graph. An alternative parametrization consists in using one dihedral angle and the twist angle τ e between the edge e and the opposite edge, see Fig. 2. The latter parametrization is closely related to the Kapovich-Millson [START_REF] Kapovich | The symplectic geometry of polygons in euclidean space[END_REF][START_REF] Bianchi | Polyhedra in loop quantum gravity[END_REF] conjugate variables (τ e , µ e ), where µ e = (j 2 ab + j 2 ac -2j ab j ac cos φ In terms of these variables, we considered the following set of equiarea tetrahedra: These data correspond to five tetrahedra that, although all have the same areas, cannot be assembled to form a 4-simplex. As a result, we also have a certain freedom in picking up how to orient them in space. To keep things as simple as possible, we choose to orient all tetrahedra aligning a normal to a face with the x axis. We further fix the rotational gauge by assuming that a normal to a second face lies in the xy plane.

φ (1) 23 = π/2 φ (2) 12 = π/3 φ (3
The resulting normals are summarized in Table 2, and it can be explicitly checked that it is not possible to find local rotations so that all pairs at a given face satisfy the orientation equations [START_REF] Perez | The Spin Foam Approach to Quantum Gravity[END_REF]. Table 2: Choice of normals for the closed twisted geometry described in the text. The gauge is fixed so that for every tetrahedron there is one normal in the x direction, and a second normal in the xy plane. The additional alignment n23 = (0, 1, 0) is a consequence of the straight angle φ23 = π/2 among the shape parameters.

ñ b b 1 2 3 4 5 1 (1, 0, 0) (0, 1, 0) (-0.

A vector geometry

Vector geometries were originally characterized in [START_REF] Barrett | Asymptotics of relativistic spin networks[END_REF], and an explicit parametrization was recently provided in [START_REF] Donà | SU(2) graph invariants, Regge actions and polytopes[END_REF] Table 3: Choice of normals for the vector geometry [START_REF] Sarno | 2-vertex Lorentzian Spin Foam Amplitudes for Dipole Transitions[END_REF] described in the text. The gauge is fixed so that all normals are pairwise antiparallel, n ab = -n ba .

A Euclidean Regge geometry

By Euclidean Regge geometry we refer to boundary data that can be embedded in R 4 to form a Euclidean flat 4-simplex. To construct such data, we start from a flat Euclidean 4-simplex, as the convex envelope of its five vertices. Using the freedom to do SO(4) transformations, we can always put a vertex in the origin of R 4 , one on the x axis, a third one in the xy plane, and a fourth one in the xyz hyperplane. This exhausts the rotational freedom, and the position of the last vertex is free. From the vertices we compute the edge vectors I e , with I = 0, 1, 2, 3 a Cartesian coordinate index. Their lengths e := || I e || := I e J e δ IJ provide ten numbers that characterize uniquely the 4-simplex up to SO(4) transformations. From the edge vectors we can also define a simple bivector to each triangle t of the 4-simplex, given by B IJ t := 2

[I e J]
e , where ee is any pair of edges in t. Then from A 2 := ||B IJ || 2 = B IJ B IJ /2, or equivalently using Heron's formula from the edge lengths, we obtain the ten areas of the 4-simplex. 7The amplitude's boundary data are areas and 3D normals in the time-gauge frame. Therefore the next step is to find the transformation that takes all tetrahedra of our 4-simplex to the time-gauge, then compute the resulting 3D normals to the faces. Only then we can identify the SO(4)-invariant geometric area with the SU(2) irrep label j. To do so, we take as reference frame the t = 0 hyperplane to which the first tetrahedron belongs (the one determined by the first four vertices fixed earlier). Its 4D normal is already the canonical one for the time-gauge, t I = (1, 0, 0, 0). Hence we can directly compute the areas and 3D normals of this tetrahedron from the edge vectors, and take them as the first set of boundary data. For each of the remaining four tetrahedra, the procedure is not as straightforward. First, we determine the 4D normal from the triple product of edges with a common vertex, e = 1, 2, 3,

N aI = IJKL J a1 K a2 L a3 || IJKL J 1 K 2 L 3 || , ( 12
)
ordered so that it is outgoing. 8 The SO(4) transformations mapping these vectors to t I can be found with a formula similar to Rodrigues's 3D rotation formula,

Λ a I J = δ I J - 1 1 + N a • t N I a N aJ + t I t J + N I a t J -(1 + 2N • t)t I N aJ , (13)
Λ a I J N J a = t I , det Λ a I J = 1.
The 4D normals define the 4D dihedral angles θ ab = arccos (N a • N b ), and provide alternative expressions for the face bivectors as well. Using the Hodge dual := (1/2) IJKL , we have in fact

B ab = ||B ab || N a ∧ N b || N a ∧ N b || = -B ba . ( 14
)
The simplicity of the bivectors is now manifest in that

N aI B IJ ab = 0 = N bI B IJ ab , (15)
and sums on the triangles in the same tetrahedron close,

b =a B ab = 0 ∀a. (16)
Using Eq. ( 13) we map each tetrahedron to the time-gauge and define the rotated bivectors in the frame of the reference tetrahedron,

• B IJ ab := Λ I aK B KL ab Λ J a L . (17)
This is simple in the time-gauge frame, i.e. t I

•

B IJ ab = 0. Therefore its magnetic part vanishes,

• B 0i ab = 0, and we can identify a 3D vector with its electric part,

j ab n k ab := k ij • B ij ab , || • B ab || = j ab . (18)
Equation ( 18) defines the remaining boundary data for the tetrahedra 2 to 5. Notice that we can also give a covariant formula for all 3D normals,

n I ab := (0, n ab ) = 2 || • B ab || t J ( • B ab ) IJ = Λ a I J N J b -N J a (N a • N b ) 1 -(N a • N b ) 2 , ( 19
)
where the case a = 1 has Λ 1 = 1 and N 1 = t. The covariant equations are actually more convenient to implement in the numerical algorithm because they use only the 4D normals, and bypass the explicit reconstruction of the bivectors. Pictorially, we can describe the above procedure as the 4-simplex being opened up by SO(4) rotations of four tetrahedra to the same R 3 frame of the first tetrahedron. The resulting 3D object was called a spike in [START_REF] Donà | SU(2) graph invariants, Regge actions and polytopes[END_REF], to which we refer the reader for explicit figures. Since the data (j ab , n ab ) defined by Eq. ( 18) or Eq. ( 19) come from a Euclidean 4-simplex, we expect that they satisfy the saddle point conditions ( 2) and (3) of the vertex amplitude asymptotic analysis. Indeed, the closure conditions (2) are satisfied since Eq. ( 16) is SO(4)-invariant, thus it applies to • B ab as well and in turn to the 3D normals. As for the orientation conditions (3), they follow from Eqs. ( 14) and [START_REF] Han | Asymptotics of Spinfoam Amplitude on Simplicial Manifold: Lorentzian Theory[END_REF]. To see this, one notices that the 3D vector (18) also coincides with the electric part of the self-dual part of • B ab , and that under Eq. ( 17) the self-dual part transforms like a rotation.

It is convenient to perform these rotations, so to have a configuration with pairwise-opposite normals (4). This configuration was called twisted spike in [START_REF] Donà | SU(2) graph invariants, Regge actions and polytopes[END_REF]. It was there shown that the required rotations have directions n 1a and angles given precisely by the 4D dihedral angle θ 1a = arccos (N 1 • N a ). In other words, the twisted spike twists the axis of the tetrahedra by an amount corresponding to the original 4D dihedral angle, and in doing so one achieves the condition of pairwise-opposite normals (4). See [START_REF] Donà | SU(2) graph invariants, Regge actions and polytopes[END_REF] again for illuminating pictures. The twisted spike has a computational advantage because all holonomies at one of the two critical points are the identity. We implemented the algorithm described above in the Mathematica code EuclideanBoundaryDataMaker in [START_REF]The code sl2cfoam is public and available at the[END_REF], which can be used to generate data corresponding to arbitrary Euclidean 4-simplices. For our numerics, we chose an equilateral 4-simplex. The main advantage is that having all areas equal, we can obtain a maximal number of data points for a fixed maximal numerical effort. The resulting normals corresponding to the twisted spike are summarized in Table 4. As a final comment it is useful to understand the relation between our "deconstruction" of the 4-simplex and the reconstruction that one does in the saddle point approximation, the critical point equations for the Euclidean spike admit the following four solutions (see Appendix B),

ñ b b 1 2 3 4 5 1 (0,
h (c) a = ± exp i 2 θ 1a n 1a σ , h (Pc) a = ± exp - i 2 θ 1a n 1a σ = h (c) a † , (20)
all with h a ∈ SU(2) ⊂ SL(2, C). These group elements are precisely the Λ a appearing in our procedure, namely one can show that

Λ a I J = 1 2 Tr σ I h (c) a σ J h (c) a † , (21)
where σ I = (1, σ), and σ I is lowered with the Euclidean metric.

A Lorentzian Regge geometry

To obtain the Lorentzian boundary data, we follow the same strategy as the previous case, but some additional care will be needed to take into account the presence of future pointing and past pointing timelike vectors. We start from a Lorentzian 4-simplex with all tetrahedra spacelike9 , and derive the corresponding 3D data in the time-gauge from it. We consider the five vertices in R 4 with metric η = diag(-+ ++), and use the SO(1,3) freedom to place them as before, a vertex in the origin, one on the x axis, a third one in the xy plane, and a fourth one in the xyz hyperplane, with the last one completely free. The tetrahedron described by the first four vertices will always be spacelike, whereas the nature of the other four tetrahedra depends on the location of the fifth vertex. We derive as before the edge vectors I e , bivectors B IJ ab , and 4D normals N I a . The norms and scalar products are now given by η, and for the epsilon symbol we take the convention 0123 = 1.

To find a convenient choice of data, notice first that a Lorentzian 4-simplex cannot be equilateral, just like a triangle in 2D. Just like in 2D, the most regular 4-simplex with all spacelike sub-simplices is isosceles, with one equilateral tetrahedron and four equal isosceles ones, so in the simplest configuration we have four areas of one value and six of a different one. To find suitable ones, we need both values to be integers or half-integers. Recall also that the numerical algorithm becomes the more costly the higher the spins involved. Therefore to optimize the numerical calculations we want a configuration in which the ratio between the two areas is maximal, while at the same time minimizing the smallest integer or half integer realization of each area value. Using the algorithm FromAreasToVericesLorentzian in [START_REF]The code sl2cfoam is public and available at the[END_REF], we scanned the space of admissible configurations, and settled on the ratio 2/5 anticipated above. 10 This is a 4-1 configuration, namely four of the timelike normals are future pointing and one past pointing, or vice versa.

Let us choose the gauge-fixed tetrahedron 1 to be past-pointing, i.e. N 1 = -t. The remaining 4 normals are all future pointing, and can be computed from the edge vectors as in [START_REF] Livine | Consistently Solving the Simplicity Constraints for Spinfoam Quantum Gravity[END_REF]. To determine the 3D normals we proceed as in the Euclidean case, with Λ a ∈ SO [START_REF] Engle | LQG vertex with finite Immirzi parameter[END_REF][START_REF] Perez | The Spin Foam Approach to Quantum Gravity[END_REF] this time. To transform the future pointing 4normals to t we select the pure boost in the plane determined by t and N a , which can be computed to be

Λ a I J = δ I J + 1 1 -N a • t N I a N aJ + t I t J + N I a t J -(1 -2N a • t)t I N aJ , a = 1, (22)
Λ a I J N J a = t I , det Λ a I J = 1.
We then have the expression [START_REF] Barrett | Asymptotic analysis of the EPRL four-simplex amplitude[END_REF] for the bivectors and Eq. ( 17) for the time-gauge bivectors as before, and define again the 3D normals via Eq. ( 18). The only difference is at the level of the covariant formula, which now reads

n I ab := (0, n ab ) = 2 || • B ab || t J ( • B ab ) IJ = -Λ a I J N J b + N J a (N a • N b ) (N a • N b ) 2 -1 , a = 1. ( 23
)
The 3D normals of the tetrahedron 1 are already in the chosen reference R 3 , and can then be directly computed from the vertices (or equally from [START_REF] Bianchi | Graviton propagator in loop quantum gravity[END_REF] but with now Λ 1 = -1, due to the fact that the normal to this tetrahedron is past-pointing). The resulting 3D object can be referred to as an antispike, because the tetrahedra boosted in the frame of the first will be "sitting inside" it, see Fig. 3 for a 1+1 example. And notice that the normals of the faces shared with the first tetrahedron will be aligned, and not antialigned. While this configuration is a good representation of a Lorentzian 4-simplex in 3D, it will not feed a critical point to the coherent amplitude [START_REF] Ruhl | The Lorentz Group and Harmonic Analysis[END_REF], because of the minus signs in front of half of the vectors, which were chosen for convenience in dealing with the vector and Euclidean geometries. 11 This can be fixed if we turn around the reference tetrahedron applying an inversion, namely the composition T P ∈ SL(2, C) of a time reversal and a parity transformation. Doing so we obtain again a picture with all tetrahedra outside, which we call the Lorentzian spike. Pictorially, it can be distinguished from a Euclidean spike because the tetrahedra spiking out are too "short" to allow the mapping to a closed Euclidean 4-simplex, see again Fig. 3 to get some intuition. For the Lorentzian spike, all 4-normals are future pointing, and the 3D normals of the first tetrahedron are each antialigned with the corresponding one of the tetrahedron sharing that face, as in the Euclidean We found an upper bound at 1/ √ 6 .41 > 2/5 = .4 for which the 4-simplex becomes degenerate, and a lower bound at 1/3 for which the isosceles tetrahedra become null. 11 The same would be true had we worked with the antipodal spinor conventions of [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF] instead of the antipodal vectors.

case. On the other hand, there is not now a set of five rotations so that Eq. ( 3) holds, because the Λ a are pure boosts for a = 1. This means that there is no twisted spike configuration, and no configuration with all critical holonomies at the identity. The 3D normals satisfy instead Eq. ( 5), namely they can be antialigned by a set of five complex rotations. The Mathematica code implementing this algorithm is called LorentzianBoundaryDataMaker in [START_REF]The code sl2cfoam is public and available at the[END_REF]. The normals obtained for the 2/5 ratio 4-simplex are summarized in Table 5 Table 5: Choice of normals for the Lorentzian Regge data. These are reconstructed from the chosen isosceles 4-simplex with area ratio 2/5, using the algorithm described in the text, and correspond to the spike configuration.

In determining the analytic asymptotic formula to be compared with this configuration, let us recall the definition of the Lorentzian dihedral angles θ L ab ≥ 0. Following [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF], the case with both timelike normals inside the same causal patch is called thick wedge, and the case with the timelike normals inside opposite causal patches is called thin wedge. With signature -+ ++, the scalar product of normals is negative in the first case, and positive in the second. Accordingly, we define thick wedge

N a • N b = -cosh θ L ab , (24)
thin wedge

N a • N b = cosh θ L ab . (25)
The 4-simplex of Table 5 is of type 4-1, and we have four thin angles θ L 1a and six thick ones θ L ab for a, b = 1. The solutions of the critical point equations for the Lorentzian spike are

h (c) a = ± exp (θ L 1a + iπ) n 1a • σ 2 , h (Pc) a = ± exp (-θ L 1a + iπ) n 1a • σ 2 = h (c) a -1 . ( 26
)
They are a combination of a boost given by the thin angle and an additional rotation by π in the same direction. The origin of this rotation is the inversion performed to have the data describing a spike, so that the 3D normals in the reference tetrahedron are the opposite of the corresponding ones in the adjacent tetrahedra, and not parallel as in the antispike configuration with the reference tetrahedron being past pointing. To recover the pure boost [START_REF] Freidel | From twistors to twisted geometries[END_REF] we remove this additional rotation, and indeed one can check that

Λ a I J = 1 2 Tr σ I ĥ(c) a σ J ĥ(c) a † , ĥ(c) a = ± exp θ L 1a n 1a • σ 2 , (27)
where σ I = (1, σ), and where σ I is lowered with the Minkowski metric.

4 Methods: Shelled sums and a useful approximation

The finiteness of the vertex amplitude [START_REF] Engle | Regularization and finiteness of the Lorentzian LQG vertices[END_REF] guarantees that the internal sums in Eq. ( 7) converge. However, the speed of the convergence depends on the boundary data considered, since for different configurations, both the booster functions and the {15j} symbol can decrease either exponentially or polynomially. This speed determines crucially the efficiency of our numerical method: If the convergence is fast we can keep the lowest order terms only, but if the convergence is slow one has to add up more and more terms, and the numerical algorithm becomes slower and slower.

To study the convergence, we introduce a homogeneous integer cutoff ∆ on the range of the internal spins, ∞ l ab =j ab → j ab +∆ l ab =j ab [START_REF] Dona | Asymptotics of lowest unitary SL(2,C) invariants on graphs[END_REF] and we compare the value of the amplitude for successive truncations. Notice that each of these truncated sums, or shelled sums to use the terminology of [START_REF] Dona | Numerical methods for EPRL spin foam transition amplitudes and Lorentzian recoupling theory[END_REF], can have a priori arbitrarily high internal spins l. It is only the difference with the boundary spins j that is being truncated, namely the number of terms being summed over in Eq. ( 6). The truncation ∆ = 0 gives the simplified EPRL model defined in [START_REF] Speziale | Boosting Wigner's nj-symbols[END_REF].

Let us estimate the number of terms in a given truncation, for the equispin case j ab = j. The l configurations compatible with Eq. ( 28) are (∆ + 1) 6 . For each l's configuration, we have to sum over the internal intertwiners k. Notice from the graphical representation (7) that each internal intertwiner always has at least one external spin. Hence increasing l has the effect of tightening the triangular inequalities for the intertwiners. Consequently, the external triangular inequalities from the j spins provide an upper bound to the number of k configurations, which is (2j) 4 when all j's are equal. To complete the estimate we have to include the 2j terms coming from the sum defining the reducible {15j} symbol used. Overall, this gives an upper bound of (∆ + 1) 6 × (2j) 5 [START_REF] Engle | The Loop-quantum-gravity vertex-amplitude[END_REF] configurations for the amplitude defined by boundary spins and intertwiners. For the coherent amplitude [START_REF] Dona | Numerical methods for EPRL spin foam transition amplitudes and Lorentzian recoupling theory[END_REF], one also has to sum over the external intertwiners i, which gives an additional (2j) 5 terms. Recall that for this configuration λ = 2j, the final upper bound estimate is

# ∆ ≤ (∆ + 1) 6 × λ 10 . ( 30
)
An exact numerical counting confirms this estimate, see the Euclidean data in Fig. 4. A numerical fit reported there shows that the neat effect of the triangular inequalities is to remove a bit more than half of the configurations.

Figure 4: Total number of terms to be summed over at a given λ, for Euclidean (blue) and Lorentzian (red) boundary data, for the simplified model (circles) and the first shelled sum (triangles). The lines are obtained through a numerical fit with the (∆ + 1) 6 λ 10 estimate from [START_REF] Engle | The Loop-quantum-gravity vertex-amplitude[END_REF]. For Euclidean data, the fitted coefficients give .38 at ∆ = 0 and .45 at ∆ = 1. These numbers give an indication of the percentage of configurations removed by the triangle inequalities.

For Lorentzian data, we cannot apply the counting leading to the estimate (29) because the spin configurations are isosceles. Nonetheless, the numerical fit shown in the figure is also with slope (∆ + 1) 6 λ 10 , and it is in qualitative agreement (a numerical fit with free exponent gives λ 9.96 for ∆ = 0 and λ 9.12 for ∆ = 1). This counting shows that for the same λ (and thus same number of data points in the asymptotic plots), the Lorentzian configurations are numerically much heavier. Generating this plot took about 24h on a standard laptop.

To give an idea of the complexity of the numerics involved, the amplitude for Euclidean data with λ = 16 and ∆ = 2 contains 10 14 configurations, each of which requires the evaluation of 10 4 integrals and five {6j} symbols. The evaluation of this data point alone took three weeks on our server.

For Lorentzian boundary data, the situation is worse, since one cannot use equal spins. Therefore for the same value of the rescaling parameter λ, higher individual spins will enter, increasing the number of terms to be cycled over. For the 2/5 ratio used as Lorentzian boundary data, we are limited at λ = 9, which corresponds to highest spins 45, beyond which we encounter numerical instabilities in the booster functions. We counted numerically all the configurations allowed by triangular inequalities, and the results are shown in Fig. 4. The most expensive data point we bought for Lorentzian data has λ = 9 and ∆ = 1, and 10 13 configurations to be evaluated.

These estimates show the high price of increasing the cutoff ∆ on the internal sums, making configurations for which the convergence is fast much more affordable. We also point out the λ 5 numerical cost coming from the use of coherent states in the boundary, which requires the additional summations over intertwiners. Hence our methods are much faster for calculations requiring orthogonal intertwiners and not coherent ones.

Convergence of the shells

Having discussed the numerical costs of the evaluations, let us now present some investigations on the convergence of the internal sums. We restricted attention to the four types of boundary data presented earlier, which will be relevant to study the asymptotic scaling. Figure 5 shows examples of convergence for γ = 1.2, and λ = 1. At this λ the convergence is quite fast for all four cases, with differences of 1% or less when increasing the truncation from ∆ = 4 to ∆ = 5. Closed twisted and vector geometries give plots undistinguishable from the one on the left panel, with numerical differences typically at 4 digits.

The convergence can, however, change quite significantly for different configurations or different values of γ. More importantly for our scopes, it can also change as we increase λ. An analysis like Fig. 5 for higher λ is unfortunately too costly and could not be performed. As a cheaper alternative, we looked at the gap ∆ 1 -∆ 0 between the first shelled sum and the simplified model. This is reported in Fig. 6. We see that for Euclidean data the gap remains bounded at 30% of the ∆ = 1 amplitude. If the same happens for the other differences ∆ 2 -∆ 1 , etc, it would mean that the speed of convergence is constant in λ. This possibility is supported by the full asymptotic study presented in the next section. In that context, we will also explain the oscillation seen in (the right panel of) Fig. 6 as the rephasing of the asymptotic formula performed by the shells.

For Lorentzian data, on the other hand, we observe a generic growth of the gap. This suggests that while the convergence was (very) fast at λ = 1, it will likely slow down as λ increases. This slower convergence for Lorentzian data as λ is increased will also be confirmed by the full asymptotic study of the next section, and suggests that the simplified model will miss important aspects of the full model for such boundary configurations.

Figure 6: Behavior of the difference Dv(j ab , n ab , ∆) = Av(j ab , n ab , ∆ = 1) -Av(j ab , n ab , ∆ = 0), for Euclidean (on the left) and Lorentzian (on the right) boundary data. The oscillating behavior in the Euclidean case is the Regge asymptotic behavior, as explained in the next Section. The growing behavior in the Lorentzian case is a clear indication that the simplified model is becoming subdominant with respect to the first shell as λ increases. This will be clarified in the next Section, where it will be shown that the simplified model misses the critical behavior and decays exponentially, while the first shell decays at power law.

Summarizing, our investigations indicate that convergence is fast at very small spins for all types of data considered, but that while it remains fast for twisted, vector and Euclidean Regge geometries, it slows down for Lorentzian Regge geometries.

Approximating the coherent states

In an effort to reduce the numerical cost of the coherent amplitude, we studied the effect of truncating the sums over the boundary intertwiners i. In fact, c i ( n f ) is roughly a Gaussian distribution in i [START_REF] Livine | A New spinfoam vertex for quantum gravity[END_REF][START_REF] Freidel | Holomorphic Factorization for a Quantum Tetrahedron[END_REF], see e.g. Fig. 7, the more accurate and the sharper for larger spins, and one can explore the impact of cutting the tails. An analytic study of the Gaussian width shows that it depends on both the spins and the normals [START_REF] Freidel | Holomorphic Factorization for a Quantum Tetrahedron[END_REF]. However the use of those formulas is not easy beyond the equilateral case. For our purposes it is more practical to estimate σ(j f , n f ) numerically, as explained in the caption of Fig. 7.

Figure 7: Left panel: The distribution of a coherent state with j ab = 10 and normals corresponding to an equilateral tetrahedron, as a function of the intertwiner label i = j12. To estimate numerically the width we interpolate the data points with a Gaussian, and use the function gsl stats sd of the GNU Scientific Library [START_REF] Contributors | GSL -GNU scientific library -GNU project -free software foundation (FSF)[END_REF]. Numerically it is then convenient to define a cutoff α excluding terms for which ci( n f ) ≤ σ(j, n)/α. This means cutting the tails at √ 1 + 2 log α σ(j, n). The horizontal line in the figure corresponds to α = 3, namely 1.79σ for which 93% of the Gaussian is covered. Right panel: The distribution of an isosceles coherent state with three spins equal to 10 and one equal to 25. For this configuration there are only 6 intertwiner spins allowed, and the spread is too broad to allow any reliable truncation. These are the tetrahedra entering the most numerically economical Lorentzian boundary data.

There is however no guarantee that truncating the tails of the Gaussians gives a valid approximation for the vertex amplitude because the intertwiner dependence of the {15j} symbol can compensate for the Gaussian damping. This happens for instance in the equilateral configurations used for Euclidean boundary data, where the coherent state distribution is roughly peaked around the middle value of the intertwiner label, 12 whereas the {15j} is peaked on extremal intertwiners 0 or 2j. In this case, it turns out that cutting the tails even just a bit introduces large errors. For the simplified model and λ around 10, one needs at least 5σ in the Gaussians to reduce the error on the full amplitude below 5%. This makes the approximation basically useless and pushed us to devise an alternative scheme. Instead of truncating all five Gaussians of the vertex amplitude, we kept two Gaussians exact and truncated only the other three (it does not matter which because of the symmetry of the problem). Using the simplified model as trial, we found that truncating three Gaussians at 1.79σ (namely covering 93%) we reduce drastically the number of terms included in the summation, see Fig. 8, while introducing a small error ≤ 3%. This reduces the required computational time to roughly 1/20, for the simplified model and spins of order 10. Right panel, the improvement on the evaluation times: We compute the ratio between the number of terms Nα=3 surviving in a truncation with α = 3 and the total number Ntot. The impact of the truncation grows roughly linearly in λ (the solid line is 13.6λ -3/2 ), and we see that for instance at λ = 20 only 5% of the terms survive.

The advantage of such a truncation becomes even more important when adding the internal sums. To make sure that the error introduced is not significantly different than in the simplified model, we performed an analog study for ∆ = 1 and λ ∼ 10, reported in Fig. 9. It confirms that α = 3 is a good compromise between cutting as many terms as possible and keeping the error small. Having so chosen the truncation α = 3, we apply it to higher spins where the evaluation times are the longest. We have no quantitative estimate of the error made, but the coherence of the plots shown in the next section reassures us that we are not deviating significantly from the few percentage points established for λ ∼ 10, and the evaluation times are very conveniently shortened.

On the other hand, a similar truncation does not work for Lorentzian boundary data because of the isosceles tetrahedra there used. For the isosceles configuration there is a smaller range of intertwiners and the Gaussians are more spread out, see Fig. 10. We found that any truncation in this case introduces large errors. Therefore for Lorentzian data we contented ourselves to truncate the only equilateral coherent intertwiner, again with α = 3. This gives an improvement of roughly 1/4 of the computational time for λ = 9. On the other hand, we can introduce the truncation already at λ = 3, something that with three Gaussians truncated as for Euclidean data was introducing a bigger error. 8 for the Lorentzian configuration, truncating only on one coherent state (the equilateral one). We can cut only a single tetrahedron due to the fact that, for low rescaling parameters, the isosceles tetrahedra do not approximate a Gaussian, as shown in Fig. 7. In this case, the fit gives 0.8λ -1/2 . Right panel: The truncation relative error in the simplified model contribution for Lorentzian boundary data, used in Section 5.4.

Summarizing, in the numerics we used the Gaussian truncation with α = 3 on three coherent intertwiners for vector and Euclidean Regge data, starting at λ = 8; and a similar truncation with α = 3 but on the unique equilateral coherent intertwiner for Lorentzian data, starting at λ = 3. This approximation enters the plots shown in Fig. 6, etc. An estimated error bar of 3% could be added but would not be visible.

Numerical analysis of the asymptotic formula

In this section, we present a selection of our numerical data. We consider the four different types of boundary geometries constructed in Sec. 3, plus a completely random configuration not satisfying even closure, to test the different possible asymptotic behaviors derived in [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF].

1. For arbitrary areas and normals, or for the closed twisted geometry: the data confirm the exponential decay;

2. For the vector geometry: The data confirm a power law decay λ -12 and no oscillations;

3. For the Euclidean Regge geometry: The data confirm a power law decay λ -12 and an oscillation with a frequency given by the Regge action;

4. For the Lorentzian Regge geometry: The data indicate a power law decay λ -12 and an oscillation with a frequency depending on γ, but more conclusive results will require additional computational power.

Generic data, open and closed twisted geometry

Open boundary data

For generic data we expect from [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF] an exponential decay, namely

A v = o(λ -K ) ∀ nonnegative integer K. (31)
To confirm this numerically, we took all spins to be equal, but the normals were randomly generated. We did not truncate the coherent states, since the Gaussian profile of Fig. 7 is lost for open configurations. 13We evaluate the amplitude in the simplified model, ∆ = 0, and the first shelled summation, ∆ = 1. The vertex amplitude is exponentially suppressed, see Fig. 11. Right panel: Closed twisted geometries data, with the normals satisfying closure but not the vector geometry conditions. The line is 3.2 × 10 -12 λ -12 , added to help the eye identifying the exponential falloff.

Closed twisted geometry

A closed configuration is enough to have a critical point for the norm of the coherent intertwiners [START_REF] Livine | A New spinfoam vertex for quantum gravity[END_REF], but not for the full amplitude (6), and [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF] predicts again an exponential falloff [START_REF] Christensen | Numerical evidence of regularized correlations in spin foam gravity[END_REF]. We computed numerically the amplitude using the configuration constructed in Sec. 3.1, and verified the exponential suppression for the simplified model and the first shelled sum, see Fig. 11. In generating these data, we did not use any truncation on the coherent states. The coherent state distributions are in fact too broadly spread for this configuration, and it is thus not possible to efficiently truncate them for the values of spins considered.

Vector geometry

For vector geometries, namely data satisfying both closure (2) and orientation (3) constraints, there is one critical point (doubly degenerate), and one finds the following power law asymptotic behavior [5]

A v = 1 λ 12 N c e iλΦc + O(λ -13 ). ( 32
)
Here N c contains the (inverse square root of the) Hessian determinant, the integration measure evaluated at the critical point, and the factors of 2 and π coming from the critical point degeneracies and the Gaussian integrations. The action at the critical point is purely imaginary, hence the phase Φ c in Eq. [START_REF] Johansson | Fast and accurate evaluation of wigner 3j, 6j, and 9j symbols using prime factorisation and multi-word integer arithmetic[END_REF]. See [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF] and our appendixes for details and explicit formulas. With the vector geometry constructed in Sec. 3.2 one finds N c = 7.86 × 10 -7 .

The value of the Φ c is irrelevant because of the global phase ambiguity discussed earlier, and we restrict attention to the absolute value of the amplitude.

For the numerical evaluation of the amplitude, we fixed γ = 1.2, and evaluated first the simplified model, ∆ = 0, and then the first shelled sum, ∆ = 1. In both cases we included a truncation on three Gaussians, for reference, on tetrahedra 2, 3 and 4. Even if the configurations are not equilateral, equal area is enough to have sufficiently peaked Gaussians.

The results are reported in Fig. 12, together with the analytic prediction. We find excellent agreement for the power law. As for the numerical coefficient, the simplified model comes out a bit short, and the first shell gives an order-one contribution approaching the predicted analytical value. Notice also the short scale nonmonotonic behavior, which is introduced by the higher order corrections to the saddle point. Presumably the first shell is already converging quite well to the analytic value [START_REF] Collet | A (simple) expression of the unitary-irreducible SL(2,C) representations as a finite sum of exponentials[END_REF], and the contributions from the second and higher shells will be small. This would mean that the behavior observed at λ = 1 in Fig. 5 (left panel) is maintained at higher λ's. Left panel: The absolute value of the vertex amplitude with a vector geometry boundary for an homogeneous rescaling of the spins in log scale. The amplitude is following the power law [START_REF] Johansson | Fast and accurate evaluation of wigner 3j, 6j, and 9j symbols using prime factorisation and multi-word integer arithmetic[END_REF]. We plot the value of the amplitude in the simplified model (∆ = 0) with empty diamonds, and the first shelled summation (∆ = 1) with full triangles. Here we used the cutoff on the coherent states discussed earlier, and the whole plot took roughly 30h. The line plotted is the analytical formula. Right panel: The absolute values rescaled by λ 12 , showing a short scale nonmonotonic behavior, coming from the higher order saddle point corrections.

Euclidean Regge geometry

For Euclidean Regge data, the integral has two distinct critical points, and the asymptotic formula given in [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF] reads

A v = (-1) χ λ 12 e iλΦc N c e iλSR + N P c e -iλS R + O(λ -13 ). ( 34
)
Here χ is an integer depending on the eventual braiding of the intertwiners; c and P c the two critical points, at which the action is purely imaginary and equal up to sign, up to a global phase Φ c ;14 N c contains the (inverse square root of the) Hessian determinant, but also a contribution from the integration measure at the critical point which is not trivial, and for convenience also the various factors of 2 and π coming from the twofold degeneracy of the critical points and the Gaussian integrations; and most importantly,

S R (j ab , n ab ) := a<b j ab θ ab ( n ab ), (35)
is the Regge action. Here θ ab are the 4D dihedral angles between tetrahedra, which can be reconstructed from the 3D normals using the spherical cosine laws.

It is always worthwhile to remind the reader that the actual Regge action is a function of edge lengths, unlike [START_REF] Martin-Dussaud | A Primer of Group Theory for Loop Quantum Gravity and Spin-foams[END_REF]. Nonetheless, the bivector reconstruction theorem [START_REF] Barrett | Relativistic spin networks and quantum gravity[END_REF] used in [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF] to derive Eq. ( 34) guarantees that the data span all Euclidean 4-simplices, and therefore Eq. ( 35) is equivalent to the Regge action. One can also give a different proof that does not rely on the bivector reconstruction theorem [START_REF] Donà | SU(2) graph invariants, Regge actions and polytopes[END_REF]. To that end, notice first that due to the rotational invariance at each node of the functions θ ab , the independent variables of ( 35) are areas and angles between the 3D normals. It is known from [START_REF] Dittrich | Area-angle variables for general relativity[END_REF] that the area-angle action ( 35) is equivalent to the Regge action based on edge lengths provided closure and shape-matching conditions are satisfied, a result that holds for a generic triangulation and not only a single 4-simplex. In the case at hand of a single 4-simplex, closure is guaranteed by the critical point condition [START_REF] Freidel | A New Spin Foam Model for 4D Gravity[END_REF]. And it can be shown explicitly [START_REF] Donà | SU(2) graph invariants, Regge actions and polytopes[END_REF] that the conditions for the existence of two distinct critical points is precisely the shape-matching conditions, in the form of consistency of spherical cosine laws. 15If N c = NPc , the asymptotics can be written as a cosine of λS R , with a λ-independent phase offset given by the argument of N c ,

A v = (-1) χ λ 12 2|N c |e iλΦc cos λS R - 1 2 arg H c + O(λ -13 ). ( 36
)
We are not aware of an analytic proof of this property, but we verified it in all cases numerically checked.

To avoid confusions between the phases Φ c and arg H c , we will refer to the first as the global phase and to the second as the cosine phase offset. For our numerical simulations, we considered the equilateral twisted spike configuration described in Sec. 3.3. Assigning the normals to the intertwiners according to the geometric picture of a 4-simplex, we have no braiding and χ = 0. For the Hessian, we use the explicit formulaq in Appendix B.1, and the Mathematica algorithm EuclideanHessian of [START_REF]The code sl2cfoam is public and available at the[END_REF] to evaluate it at the chosen configuration. The result is real and identical at both critical points, and fixing for example γ = 1.2, we find

N c = N P c = 2.65 × 10 -7 . (37)
Furthermore, using our phase convention for the coherent amplitude (which is not the same as the one in [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF], see Appendix A) and the conventional phase convention for the coherent states [START_REF] Perelomov | Generalized coherent states and their applications[END_REF], we numerically found that Φ c = S R . Removing this phase (redefining for instance the phase of the coherent states), the analytic prediction is

A v = 5.13 × 10 -7 λ 12 cos 5λ arccos - 1 4 . (38)
We recall that for these boundary data λ = 2j. Notice that there is no phase offset in the argument of the cosine, unlike for SU(2) asymptotics [START_REF] Ponzano | Semiclassical limit of Racah coefficients[END_REF][START_REF] Esterlis | Poisson brackets, and singular differential forms[END_REF][START_REF] Donà | SU(2) graph invariants, Regge actions and polytopes[END_REF]. This fact will play a role below.

In the numerics we used the approximation truncating three of the Gaussians, for reference the tetrahedra 2, 3 and 4. We evaluated first the simplified model, then added the internal sums, with truncations at ∆ = 1 and ∆ = 2. This increases progressively the cost of the numerics. We were able to reach λ = 40 for the simplified model, 20 for the first shell, and 18 for the second shell. We did not push our analysis to a third shell, because we would have been forced to stop at spins too small to test the asymptotic behavior. The data obtained are complex, but since we identified the global phase, we can remove it a posteriori and obtain real data points.

The first result we found is that the numerical data are in excellent agreement with the λ -12 power law, see Fig. 13. The agreement is achieved already at small spins, and holds for both the simplified model and the first shelled sum. We plot also the difference ∆ 1 -∆ 0 to show that it follows the same power law: not only does the simplified model see the critical behavior, but so do (at least some of) the individual higher l configurations.

Figure 13: Numerical data versus the analytic asymptotics [START_REF] Barrett | Asymptotics of relativistic spin networks[END_REF], for Euclidean Regge data, γ = 1.2. Log-log plots, absolute values of the amplitude, versus the analytic scaling 5.13 × 10 -7 λ -12 . Left panel: Simplified model and first shelled sum. Both scale according to the predicted power law, and feature oscillations. Right panel: The simplified model and the difference between the first shelled sum and the simplified model; namely, all contributions with one nonminimal l. This difference also scales with the predicted power law and features oscillations, showing that individual higher l contributions see the critical behavior.

We found the same λ -12 scaling also for the second shell ∆ 2 -∆ 1 -∆ 0 , with includes individual terms with two l = j + 1 or one l = j + 2. We deduce from these numerical results that the simplified model captures the right scaling, and so do individual shells. This means that the convergence of the internal sums for higher λ is still qualitatively similar to the extended study that was possible to make for λ = 1 (see Fig. 5), and also that the the analytic behavior (34) must be the result of adding up all the relevant contributions of different l's. These are of order one initially, but we expect them to suitably decrease as guaranteed by the convergence.

Next, let us zoom in on the palatable oscillations already visible in Fig. 13. We rescale the data and asymptotic formula by λ 12 and report the result in Fig. 14, the highlight plot of our paper. It shows a beautiful agreement between the numerical evaluation of the vertex and the analytic asymptotic formula in particular, a confirmation of the frequency of the oscillations determined by the Regge action, its most important feature. 16Figure 14: Rescaled numerical data versus the analytic asymptotics [START_REF] Barrett | Asymptotics of relativistic spin networks[END_REF], for Euclidean Regge data and γ = 1.2, for the simplified model and internal summations with cutoffs ∆ =1 and 2. The simplified model alone already sees the Regge oscillations of the analytic formula. The two shelled sums maintain the same frequency, but give contributions of order one to the magnitude and the phase offset.

In more details, we observe the following situation. First, the simplified model captures already the right frequency of oscillations. What is missing is the precise magnitude and cosine phase offset. The first and second shells have the same frequency of oscillations, and different magnitude and cosine phase offset, which lead to a better qualitative match with the asymptotic formula. A more precise quantitative improvement will likely require pushing to higher spins where the asymptotic formula becomes more accurate. As for the crude numerical values, we see that at λ ∼ 20 (which means spins j ∼ 10, recall half-integers are also being computed here) the simplified model and shelled sums agree with one another and with the analytic formula with a 10% error.

From these data we deduce that individual l's see the same critical behavior for Euclidean Regge data, including the simplified model, and that the asymptotic formula [START_REF] Rashid | Boost Matrix Elements Of The Homogeneous Lorentz Group[END_REF] is the result of summing all the relevant contributions. However the power law and frequency of oscillations are well described by the simplified model, with the internal sums only contributing to the overall magnitude and phase offset. Given the much faster evaluation times of the simplified model, this is an interesting property, whose origin will be explained in the next section.

The data presented in the figures of this section have γ = 1.2. For double checking, we performed further numerical investigations with the different value γ = 0.1, leading to similar plots. In particular, we found the same frequency of oscillations as predicted by Eq. [START_REF] Rashid | Boost Matrix Elements Of The Homogeneous Lorentz Group[END_REF].

We conclude from the numerics that the power law and frequency of the asymptotic formula can be confirmed already at spins of order 10, to a 10% error. A better match requires changes in the amplitude and the phase offset, for which one needs more shells and/or higher spins.

Lorentzian Regge geometry

For Lorentzian Regge data, we have two distinct critical points, and the oscillating power law decay [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF] A v = (-1) χ+M λ 12 e iλΦc N c e iλSc + N P c e -iλSc + O(λ -13 ).

(

) 39
The crucial difference is that this time the action at the critical point depends linearly on γ,

S c = γ a<b j ab θ L ab ( n ab ). ( 40
)
The Lorentzian 4D dihedral angles θ L ab are defined from the Lorentzian spherical cosine laws in terms of the 3D normals, or from the 4D normals as recalled in Sec. 3.4 above. They are boosts between the timelike normals if the latter are both future or past pointing, and in the mixed case (called the thin wedge configuration) one has to add a shift by a factor iπ. These shifts result in the additional (-1) M phase, where

M = thin wedges j ab . ( 41
)
We remark that Eq. ( 40) is a Lorentzian Regge action with areas are given by γj, which is in agreement with the area spectrum of loop quantum gravity in the large spin limit. As for the shape-matching and equivalence of Eq. ( 40) to the Regge action based on lengths, the same considerations of the Euclidean case apply.

For our numerical simulations, we considered the isosceles configuration described in Sec. 3.4. Assigning the normals to the intertwiners according to the geometric picture of a 4-simplex, we have no braiding and χ = 0. Further, we have four thin wedges all associated with j ab = 5λ, thus M = 20λ with integer λ and this phase drops out as well. For the Hessian, we use the explicit formulas in Appendix B.1, and the Mathematica algorithm LorentzianHessian of [START_REF]The code sl2cfoam is public and available at the[END_REF] to evaluate it at the chosen configuration. Taking γ = 0.1, we find

N c = N P c = 1.76 × 10 -13 + i1.87 × 10 -14 , (42)
thus Eq. ( 39) is a cosine up to a global phase. Reabsorbing the global phase in the definition of the coherent states and computing the corresponding Regge action, we obtain the analytic asymptotic formula

A v = 3.53 × 10 -13 λ 12 cos (0.01λ + 0.106) . ( 43
)
In running the numerical code, we used the Gaussian truncation on only one coherent intertwiner, the equispin gauge-fixed tetrahedron 1. The data obtained are complex, and their λ-dependent phase this time is not simply the Regge action. We spent a considerable amount of time going through the conventions and technical details that determine this anyway irrelevant global phase, for the sake of having real data points. In the end the simplest solution was to determine the λ-dependent phase through a numerical fit, then remove it from the numerical data. A more brutal removal by taking absolute values would have deprived us of the numerical sign, which is important in comparing data and analytic formulas. 17 The code works fine also for Lorentzian boundary data, and we were able to evaluate the amplitude for the simplified model and first shell. The numerical results we were able to obtain are however less satisfying than for the previous configurations, in terms of comparing with the analytic asymptotic formula. The main problem is the numerical instability in the booster functions, which limits us to λ = 9, before the j = 50 instability of the boosters is reached. Furthermore, the numerical cost of the calculation limits us to ∆ = 1. We can, on the other hand, take advantage of the expected linear dependence of the oscillations on γ, and evaluate numerically the vertex amplitude at different values of γ, to try to establish different aspects of the asymptotics. We considered three different values of γ, chosen to have qualitatively different behaviors in the accessible λ ∈ [START_REF] Engle | LQG vertex with finite Immirzi parameter[END_REF][START_REF] Ponzano | Semiclassical limit of Racah coefficients[END_REF] range: a rather flat plateau with no oscillations, and visible oscillations but with a different frequency. The resulting asymptotic formulas needed to match the numerics in the three cases are reported in Table 6. 17 We shamelessly admit that we actually failed to reproduce the observed numerical global phase from our analytic calculations -Ref. [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF] does not compute it, given its irrelevance. Appendix A collects all explicit formulas and all numerous places where a global phase arises, as well as a careful comparison between the conventions here used based on [START_REF] Varshalovich | Quantum Theory of Angular Momentum: Irreducible Tensors, Spherical Harmonics, Vector Coupling Coefficients, 3nj Symbols[END_REF][START_REF] Perelomov | Generalized coherent states and their applications[END_REF][START_REF] Ruhl | The Lorentz Group and Harmonic Analysis[END_REF], and the ones of [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF]. Albeit an overlooked phase lurking there and open to inspection, we suspect that the origin of the problem lies in a mismatch in the conventions for the orientation of the links of the reducible {15j} between the numerical code and our analytic description of the amplitude. Not having enough points for a significative log-log plot, we present directly the data rescaled by the predicted λ 12 falloff, versus the rescaled asymptotic formula. The first plot shown in Fig. 15 is for γ = 0.1. The effective frequency of oscillations is low, giving a roughly constant behavior in the accessible range. The first consideration to be deduced from the data is that the simplified model appears to fall off significantly faster than the analytic power law decay of the EPRL model. We interpret this to mean that the simplified model does not have a critical behavior for Lorentzian Regge geometries. The data points of the first shell on the other hand roughly stabilize in parallel to the asymptotic formula. This indicates that the first shell is already enough to capture the critical behavior, even though with its actual magnitude being too small, we expect nearby shells to still give order-one contributions, and thus a slow convergence. Figure 15: Rescaled numerical data for Lorentzian Regge geometry, versus the analytic asymptotics [START_REF] Freidel | Holomorphic Factorization for a Quantum Tetrahedron[END_REF], for γ = 0.1. The bottleneck is the numerical instabilities of the booster functions at spins of order 50, which even for the optimal isosceles configuration we identified, occur already at λ = 10. With the caveat of the few points obtained, the simplified model is decaying faster than the first shell, indicating that it lacks the critical behavior, and will likely decay exponentially. The first shell, on the other hand, is well in line with the power law decay. Its magnitude falls short of the asymptotic formula by a factor of order one: this could be completed by higher order shells, accordingly with the previous indications that convergence is slow for Lorentzian Regge geometries.

Next, we increase the expected frequency looking at larger values of γ. Taking γ = 7 and γ = 8 allows us to fit roughly one full period of analytic oscillations within the accessible range, and with different frequencies, to see whether the data offer a compatible nonmonotonic behavior. The results are plotted in Fig. 16. The behavior of the first shell is again compatible with the power law scaling, and we clearly observe nonmonotonic behavior and dependence on γ. This is nice evidence in support of the numerical evaluation of the vertex amplitude and of the validity of the asymptotic formula. The data points at disposal are however too few and at too small λ to draw any conclusions about the frequency of oscillations like we could do for the Euclidean Regge geometry. As a somewhat optimistic remark, as we increase γ from 7 to 8 the (interpolated) analytic zero with positive slope moves to the left, and so does the (interpolated) zero of the numerical data. Figure 16: Rescaled numerical data versus the analytic asymptotics [START_REF] Freidel | Holomorphic Factorization for a Quantum Tetrahedron[END_REF], for Lorentzian Regge data and γ = 7 and 8. These values are chosen so that the analytic formula has a visible oscillation within the range numerically accessible. The faster frequency of the γ = 8 case can be inferred comparing the eyesight interpolation of the positive slope through the x axis. The data points for the first shell are well in line with the analytic power law, and show a nonmonotonic behavior qualitatively comparable with the predicted oscillations. Higher spins are needed for a better match, as well as increasing the number of shells being summed over.

In summary, the simplified model does not see the critical behavior for Lorentzian Regge data, but the first shell already does, and presumably all other shells. The power law decay observed confirms the asymptotic formula. There is a γ-dependent nonmonotonic behavior, qualitatively compatible with the oscillations predicted by the asymptotic formula. The convergence of the internal sums is slower than with Euclidean Regge data, and higher shells/higher spins are needed for a more quantitative confirmation of the asymptotic formula.

Outlook

Geometric meaning of the internal sums

The formula [START_REF] Dona | Numerical methods for EPRL spin foam transition amplitudes and Lorentzian recoupling theory[END_REF] and more generally the method of [START_REF] Speziale | Boosting Wigner's nj-symbols[END_REF], introduce a factorization of Lorentzian spin foam amplitudes so that only SU(2) objects appear at the vertices, and all Lorentzian properties are in the booster functions localized at the edges. Since the vertices are dual to 4D structures and the edges to 3D structures, this raises the question of how the 4D Lorentzian structures are reconstructed by the internal sums. A second question concerns the simplified model, which we have shown to capture the EPRL asymptotics for vector and Euclidean Regge geometries. Since the simplified model has the advantage of much faster evaluation times, it is useful to understand why this is the case. These two questions are closely related, and can be answered by inspecting in more detail the role of the sums over the internal spins l ab .

Let us look again at the decomposition [START_REF] Dona | Numerical methods for EPRL spin foam transition amplitudes and Lorentzian recoupling theory[END_REF], and consider the simplified model, for which the internal sums are removed. 18 The vertex amplitude reduces to the {15j} symbol evaluated at the boundary spins, and the booster functions at their minimal configurations:

A EPRLs v (j ab , n ab ) = ka,ia {15j} 5 a=2 d ka B γ 4 (j ab , j ab ; i a , k a ) 5 a=1 c ia ( n ab ). (44)
The internal intertwiners are still being summed over, but in the large spin limit the booster functions are Gaussians peaked on equal intertwiners, and to lowest order [START_REF] Puchta | Asymptotic of Lorentzian Polyhedra Propagator[END_REF][START_REF] Speziale | Boosting Wigner's nj-symbols[END_REF][START_REF] Sarno | 2-vertex Lorentzian Spin Foam Amplitudes for Dipole Transitions[END_REF]]

B γ 4 (j ab , j ab ; i a , k a ) b(γ, j ab ) λ 3/2 d ia δ iaka . (45)
The simplified vertex amplitude [START_REF] Ruhl | The Lorentz Group and Harmonic Analysis[END_REF] is thus in first approximation proportional to a single coherent {15j} symbol, which we recall is the vertex amplitude for SU(2) BF theory:

A EPRLs v (j ab , n ab ) 1 λ 6 5 a=2 b(γ, j ab ) A SU(2) v (j ab , n ab ), A SU(2) v (j ab , n ab ) := ia {15j} 5 a=1 c ia ( n ab ). ( 46
)
This approximation shows that the simplified model has the same critical points of SU(2) BF theory, and λ -6 times its scaling. The SU(2) coherent vertex amplitude has no critical behavior for Lorentzian Regge data and thus decays exponentially; whereas for Euclidean Regge data it admits two distinct critical points and oscillations with exactly the same frequency [START_REF] Martin-Dussaud | A Primer of Group Theory for Loop Quantum Gravity and Spin-foams[END_REF] of the EPRL model [START_REF] Barrett | Quantum gravity asymptotics from the SU(2) 15j symbol[END_REF],

A SU(2) v (j ab , n ab ) = 1 λ 6 1 2 6 π 2 e iλΨc | det H SU(2) c | 1/2 cos λS R - 1 2 arg H SU(2) c + O(λ -7 ). (47)
Here H SU (2) c is the Hessian determinant at the critical point, and Ψ c a global phase determined by the action at the critical point, which is purely imaginary, and which depends in turn on the gauge choice and the phase of the coherent states exactly as for the EPRL model.

From this analysis we understand why the simplified model has no critical behavior for Lorentzian Regge data, and the same critical behavior of SU(2) BF for vector and Euclidean Regge data. In particular, for Euclidean Regge data, inserting Eq. ( 47) into Eq. ( 46) we find that the simplified model has the same λ -12 power law decay and the same γ-independent frequency of oscillations S R of the full EPRL model. This explains why the simplified model captures the right scaling and frequency of oscillations of the EPRL asymptotics [START_REF] Rashid | Boost Matrix Elements Of The Homogeneous Lorentz Group[END_REF], as was shown in Fig. 14.

It also explains why it does not capture the right cosine phase offset: The simplified model sees the SU(2) phase offset - 1 2 arg H SU(2) c = -0.324 for the equilateral configuration, which is absent in the asymptotic formula for the EPRL model, since the Hessian is in that case real. Consider then a modified asymptotic formula

A s v (Euclidean data) = 5.13 × 10 -7 λ 12 cos 5λ arccos - 1 4 -0.324 , (48)
where we added by hand the offset of the SU(2) BF amplitude. This formula correctly matches the cosine phase offset of the simplified model, see Fig. 17. As a consequence, the internal sums pile up to the SU(2) Hessian to give the right magnitude of the EPRL amplitude, but also contribute to the phase offset of the cosine, creating an interference pattern that changes Eq. ( 47) to Eq. [START_REF] Rashid | Boost Matrix Elements Of The Homogeneous Lorentz Group[END_REF]. with the SU(2) offset [START_REF] Anderson | Recursion and symmetry relations for the clebsch-gordan coefficients of the homogeneous lorentz group[END_REF]. The improved agreement shows that the right asymptotics of the simplified model is to first approximation the same of SU(2) BF theory, and implies that the role of the internal sums in this case is to contribute to the re-phasing of the offset to match the EPRL asymptotics -as well as contributing to the magnitude. Right panel: Evauation of the SU(2) BF asymptotics for an equilateral Euclidean 4-simplex, and comparison with the analytic formula (47) (It updates the plot of [START_REF] Donà | SU(2) graph invariants, Regge actions and polytopes[END_REF] which stopped at λ = 50 -or λ = 25 in that notation). It shows how while the frequency is captured early on, the magnitude requires higher spins. At λ = 79 the data for the {6j} required weighted 230 gigabytes and saturated our server's capacity.

We take this opportunity to provide in the same figure an updated plot of the SU(2) BF asymptotics with respect to the one presented in [START_REF] Donà | SU(2) graph invariants, Regge actions and polytopes[END_REF]. It allows us to highlight that the frequency of the asymptotic formula is matched early on by the exact evaluation, but the magnitude only at higher spins. We expect a similar situation for the EPRL model.

Let us now consider Lorentzian boundary data. As already stated, the SU(2) vertex amplitude is exponentially suppressed for such data and thus also the simplified model should be, in agreement with the indication of exponential decay in Fig. 15. Only the l-shells can see the critical behavior. If we look at the factorization (6), we need both booster functions and the {15j} symbol to individually have critical behavior. The {15j} has critical behavior only for Euclidean Regge configurations, or vector geometries. The booster functions must then admit a critical behavior precisely at those configurations for which the j's correspond to a spacelike Lorentzian 4-simplex, and the l's to a Euclidean one, or a vector geometry. A saddle point approximation of the booster function will appear in [START_REF] Dona | Asymptotics of SL(2,C) coherent invariant tensors[END_REF], and shows that critical behavior appears when the two sets of spins and coherent intertwiners correspond to tetrahedra that can be boosted into one another. The action of the boosts is defined embedding each 3D normal as the electric part of a γ-simple bivector. Hence there are a priori an infinite number of admissible l's contributing to the critical behavior, characterized by the existence of a certain classical boost mapping the j tetrahedra to the l tetrahedra, and then a critical point of the {15j} symbol at that configuration. There are a priori infinitely many Euclidean 4-simplices for which this construction works, and we had fun exploring the ones with smallest l's in a reconstruction algorithm depicted in Fig. 18.

According to this argument, all shells should contribute to the asymptotic formula, and Eq. ( 39) should only be reproduced accurately when sufficiently many shells are summed over. One can reasonably hope that a limited number of shells will suffice to most applications, and in the luckiest cases, the single first shell and its qualitative matching to the asymptotic behavior. For situations where one needs to sum many shells, the method [START_REF] Speziale | Boosting Wigner's nj-symbols[END_REF] here used becomes less efficient and may have to leave way for alternative approaches, like for instance adaptive Monte Carlo. 6.2 Numerical improvements

Improved server

A few of the numerics presented in this paper can be performed on a simple laptop, like computing the simplified model at spins below 10, our counting of the number of configurations compatible with the triangle inequalities. For the rest, we used a server with 32 cores at the CPT. While this enormously improves with respect to a single laptop, it is effectively very small with respect to the computing power used in other branches of theoretical physics. The CPT is in the process of acquiring a server with 128 cores. When it will come on line, it will cut our evaluation times by a factor of 4, and push forward current computing and RAM limits. To fully take advantage of this additional power however, one will need to improve the part of the code that computes the booster functions, to eliminate the numerical instabilities plaguing it at spins of order 50. Once this is done, a more powerful server will no doubt be able to fully test the asymptotics with Lorentzian boundary data, and show at which spins it is reached, and for which cutoff ∆, answering the questions left open here.

A more powerful server will also allow numerical evaluations for extended triangulations. In particular for calculations not requiring the use of coherent states, computational times may be reasonable for a triangulation with a handful or more of 4-simplices. Beyond that, one may likely need simplification schemes, like for instance the minisuperspace models used for a numerical study of renormalization in the Euclidean EPRL [START_REF] Bahr | Hypercuboidal renormalization in spin foam quantum gravity[END_REF], or approximations. To that end, let us briefly comment on one approximation and two variations that could be considered in future work.

Cutting the boosters

Cutting the coherent boundary states is the only approximation used in some of the calculations of this paper. A second approximation that should be considered is to cut the booster functions, as it could turn out to be more efficient. The boosters are decreasing functions of the l's, typically power law but in many cases exponentially decreasing. If one could identify a priori which ones, then one could cut the l sums. The numerical investigations we performed do not suggest any easy way in which this can be done. We found for example cases in which most of the boosters within a shell are one to 3 orders of magnitude smaller than the few dominant ones, yet the total sum of the subdominant ones is of the same order as the sum of the dominant ones, thus forestalling any truncation. Some light over this possibility could be shed by analytic studies of the boosters' asymptotic behavior [START_REF] Puchta | Asymptotic of Lorentzian Polyhedra Propagator[END_REF][START_REF] Dona | Asymptotics of SL(2,C) coherent invariant tensors[END_REF].

Alternative expressions for the boosters

Our code computes the booster functions using numerical integration with the trapezoidal method, and the expression of the boost matrix elements d as finite sums derived in [START_REF] Collet | A (simple) expression of the unitary-irreducible SL(2,C) representations as a finite sum of exponentials[END_REF]. This expression contains ratios of Euler's gamma functions with both real and complex arguments. We spent a significative amount of time optimizing this part of the code and running tests and debugs. The instabilities we found arising at spins of order 50 are the results of very small numbers constructed from ratios of very large numbers, and it is the part of our code that we should, and we plan to, improve next.

On the other hand, there are alternative expressions that could be considered in future work. One is the alternative formula derived in [START_REF] Speziale | Boosting Wigner's nj-symbols[END_REF], where the r integrals are solved analytically but at the price of introducing a new integral over a virtual ρ label. The ρ integrand has a nice Gaussian-like behavior, 19 , which can be exploited in numerical evaluations.

Another possibility is to implement the recursion relations for the SL(2, C) Clebsch-Gordan coefficients derived in [START_REF] Anderson | Recursion and symmetry relations for the clebsch-gordan coefficients of the homogeneous lorentz group[END_REF]. For a simpler model with three-valent boosters and no intertwiners, this approach is extremely powerful, and one can push the internal sums to order ∆ ∼ 100 [START_REF] Dona | Searching for classical geometries in spin foam amplitudes: a numerical method[END_REF].

Alternative form of the coherent states

The coherent state coefficients c i ( n ab ) are a fast part of the numerical code. However, since an enormous number of these coefficients is needed, it could be interesting to explore alternative definitions that could speed up this procedure. One possibility is the cross ratio coherent states [START_REF] Freidel | Holomorphic Factorization for a Quantum Tetrahedron[END_REF], and the analytic expressions studied in [START_REF] Bonzom | Generating Functions for Coherent Intertwiners[END_REF]. This basis has more economical features and proved advantageous in numerical evaluations of the coherent state volume operator defined in [START_REF] Bianchi | Polyhedra in loop quantum gravity[END_REF]. We hope to be able to investigate its implementation in future work.

As we estimated in Sec. 4, the use of coherent states can require up to an additional λ 5 terms, which has a strong impact on the numerical cost of the amplitude. Coherent states are necessary for the asymptotic formula, but many spin foam calculations can be done in the more economical orthonormal basis of intertwiners. This can be an important trade-off when moving to consider amplitudes on many 4-simplices. We have seen that the number of terms for a single equilateral 4-simplex with coherent states scales like λ 10 . This is the same number of terms needed to compute, with orthonormal intertwiners, the equispin amplitude with three 4-simplices glued together and one internal sum. This ∆ 3 amplitude is relevant for a discussion of the flatness problem [START_REF] Bayle | Short-scale Emergence of Classical Geometry[END_REF][START_REF] Hellmann | Holonomy spin foam models: Asymptotic geometry of the partition function[END_REF], and will be the object of future work. More in general, avoiding coherent states and increasing the number of cores, we expect to be able to make calculations on a few vertices. Going beyond that will probably require some deeper optimization of the method, or reduction to minisuperspace models.

Conclusions

We presented for the first time a numerical evaluation of the Lorentzian EPRL vertex amplitude. Obtaining this result proved to be a challenging task, and the work presented in this paper took more than two years to complete. Preliminary, more qualitative results were presented at the conference Loops '17. Making them quantitatively precise was complicated by the many technical details to be worked out. The challenges were twofold. The main one was clearly to setup the numerical code, which includes evaluating the booster functions, the {15j} symbols, the coherent states, and assembling them via hash tables. All making sure that there were no numerical instabilities or inconsistent conventions. A bit more surprisingly, we also faced some analytical challenges, concerning the way the boundary data were to be constructed and fed to the computer, and computing the overall phase and Hessian of the asymptotic formula. While these quantities were well defined in principle, precise explicit formulas and reconstruction algorithms were missing. We devoted part of the paper to explaining how this can be done explicitly, and the online repository [START_REF]The code sl2cfoam is public and available at the[END_REF] contains explicit Mathematica codes to answer these questions. We hope that our paper and this repository can provide a reference point for future work in this direction, so that one does not have to spend time again on this background material.

On the analytic side, we provided useful formulas and explicit configurations for the boundary data to be used in the amplitude. We shed light on the geometric meaning of the internal sums, explaining in particular why the simplified model asymptotes to the Regge action for Euclidean configurations but not for Lorentzian ones.

On the numerical side, our main conclusion is that precise numerical evaluations of the complete Lorentzian EPRL model are currently possible, and it is worthwhile to invest computational capabilities to address explicit physical questions. The method of [START_REF] Speziale | Boosting Wigner's nj-symbols[END_REF] and its implementation developed in [START_REF] Dona | Numerical methods for EPRL spin foam transition amplitudes and Lorentzian recoupling theory[END_REF] have been successfully tested. Including only a few terms of the internal shells appears sufficient to most applications, avoiding the need to develop a full Monte Carlo approach of the amplitude. The codes for the geometric reconstruction of the data and the numerical evaluation of the vertex amplitudes are publicly available at [START_REF]The code sl2cfoam is public and available at the[END_REF]. More specifically to the asymptotic large spin behavior of the vertex amplitude, our conclusions confirm the analytic formula of [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF] already at spins of order 10 for vector and Euclidean Regge data. Higher spins are required for a good match in the case of Lorentzian Regge data, but the power law decays and γ-dependence of the oscillations are already evident at spins of order 10.

The computational resources it takes to compute the vertex amplitude in the coherent basis are considerable but at the same time affordable. A natural and indeed crucial question is how to extend the numerical calculations to many vertices. We performed computations using 32 parallel cores, but let us not forget that in other areas of physics, computations are done with hundreds or thousands of cores. What will be computed in the future depends on the resources available, as well as on numerical improvements and optimizations that can still be done. Finally, we would like to stress that while our methods and results were motivated by the EPRL model, they can be applied or adapted to other spin foam models based on or related to SL(2, C) [START_REF] Kaminski | Spin-Foams for All Loop Quantum Gravity[END_REF][START_REF] Hellmann | Holonomy spin foam models: Asymptotic geometry of the partition function[END_REF][START_REF] Haggard | SL(2,C) ChernSimons theory, a non-planar graph operator, and 4D quantum gravity with a cosmological constant: Semiclassical geometry[END_REF][START_REF] Engle | Lorentzian proper vertex amplitude: Classical analysis and quantum derivation[END_REF][START_REF] Engle | Lorentzian proper vertex amplitude: Asymptotics[END_REF][START_REF] Kaminski | Asymptotic analysis of the EPRL model with timelike tetrahedra[END_REF][START_REF] Liu | Asymptotic analysis of spin foam amplitude with timelike triangles[END_REF], their physical applications, e.g., [START_REF] Christodoulou | Planck star tunneling time: An astrophysically relevant observable from background-free quantum gravity[END_REF], to applications involving many vertices through intertwiner renormalization as in [START_REF] Dittrich | Decorated tensor network renormalization for lattice gauge theories and spin foam models[END_REF] or minisuperspace renormalization as in [START_REF] Bahr | Hypercuboidal renormalization in spin foam quantum gravity[END_REF], as well as to calculations of tensorial and graph invariants regardless of their quantum gravity applications.

n → -n ⇔ (Θ, Φ) → (π -Θ, π + Φ) ⇔ ζ → ζ AP = - 1 ζ . (A.9) We observe that |j, ζ] = e 2ij(π+Φ) |j, ζ AP , |j, ζ AP ] = e 2ijΦ |j, ζ , (A.10)
so the antipodal map has the effect of changing family without preserving the Hopf section because of the additional phase. This will play a role below.

Coherent intertwiners

Group averaging the tensor product of coherent states defines the coherent intertwiners [START_REF] Livine | A New spinfoam vertex for quantum gravity[END_REF],

SU

dg

⊗ f D (j f ) (g)|j f , ζ f = i12 d i12 c i12 ( n f ) |j 1 , . . . j 4 , i 12 , (A.11)
where i 12 is the standard orthonormal intertwiner basis, here taken in the 12 recoupling channel. The coherent intertwiner coefficient c i12 ( n f ) is given by Eq. ( 9) in terms of Wigner's four-legged symbol

j f m f (i12) = j 1 j 2 j 3 j 4 m 1 m 2 m 3 m 4 (i12) = m12 (-1) j12-m12 j 1 j 2 j 12 m 1 m 2 m 12 j 12 j 3 j -m 12 m 3 m .
In Eq. ( 9) we used n to label the argument of the Wigner matrices, instead of Eq.(A.2). This was done so to reduce at most the technical details needed to read the main text. This coherent intertwiner is directly applicable to the tetrahedron labeled 5 in Eq. ( 7). In the other cases, one or more kets are replaced by bras, and simultaneously we put an antipodal map in the normal labeling the coherent states, like in Eq. ( 1).

Where we have a bra, the intertwiner needs to be multiplied by the epsilon matrix (j) to take into account the link's orientation. For a tetrahedron with only one bra, say in the first position of the graphic symbol, the resulting coefficient is

c (1) i12 ( n f ) := j 1 j 2 j 3 j 4 -m 1 m 2 m 3 m 4 (i12) (-1) j1-m1 D (j1) m1,j1 (-n 1 ) 4 f =2 D (j f ) m f ,j f ( n f ).
(A.12)

For the numerical code, it was convenient to rewrite this expression in terms of the intertwiner with all four legs outgoing. This can be achieved by recalling that

D (j) m,n (g) = (-1) m-n D (j)
-m,-n (g), (A.13) which applied to Eq. (A.12) leads to

c (1) i12 ( n f ) = j 1 j 2 j 3 j 4 m 1 m 2 m 3 m 4 (j12) D (j1) m1,-j1 (-n) 4 f =2 D (j f ) m f ,j f ( n f ). (A.14)
These are the expressions used in sl2cfoam.

We introduce next the antilinear map J , the infinite-dimensional analog of the complex structure J for spinors, given by

J F (z A ) := AF (z A ), (A.35a) AF (z A ) := ρ 2 + k 2 π CP 1 dµ(w A )[w|z -k-1-iρ [w|z k-1-iρ F (w A ). (A.35b) It satisfies J : (ρ, k) → (ρ, k), (J F, J F ) = (F , F ), J 2 = (-1) 2k . (A.36)
Using this map we can define the bilinear pairing

(F, F ) := (J F, F ) = (-1) 2k (F , F ), (F, F ) = (F , J F ). (A.37)
This is the infinite-dimensional analog of the spinorial bilinear (A.18).

Group elements and vertex amplitude

The SL(2, C) group action in the homogeneous representation is given by matrix multiplication on the spinor argument, h F (z A ) = F (h T z A ). The matrix elements can be equally written as the Hermitian or the bilinear products,

D (ρ,k) jmln (h) := ρ, k; j, m|h|ρ, k; l, n = F (ρ,k) jm , h F (ρ,k) ln = h F (ρ,k) ln , J F (ρ,k) jm . (A.38)
For the EPRL model, we are only interested in the "minimal" weights

ρ = γk, j = l = k, (A.39)
and for group elements in the form h = h -1 a h b . Using the invariance of the integration measure, we can write these group elements as

D (γj,j) jmjn (h -1 a h b ) = CP 1 dµ(z A )F (γj,j) jm (h T a z A )F (γj,j) ln (h T b z A ). (A.40)
The vertex amplitude in the orthonormal basis is simply the tensor associated to the complete five-valent graph,

A v (j ab , m ab ) = Because of its SU(2) invariance it can be contracted for free with the four-legged Wigner's symbol to define the amplitude in the intertwiner basis A v (j ab , i a ).

A.3 Definition of the SL(2, C) coherent amplitude

Since the maximal isotropy group of SL(2, C) is the same as of SU(2), coherent states are a direct embedding of the SU(2) ones into the unitary irreps used. For the EPRL irreps in the homogeneous representation (A.33), we define the (lowest-weight) family of coherent states

F (γj,j) jζ (z A ) := m F (γj,j) jm (z A ) j, m|j, ζ = d j π ||z|| 2(iγj-1-j) z|ζ 2j , (A.42)
where the last step follows from the factorization properties of coherent states. 20 It is useful to define also a dual basis associated with the highest weight (A.6),

F (γj,j) jJζ (z A ) := m F (γj,j) jm (z A ) j, m|j, ζ] = d j π ||z|| 2(iγj-1-j) z|ζ] 2j . (A.43)
This is related to the antipodal map by the same phase shift (A.10) of the SU(2) case, namely

F (γj,j) jJζ (z A ) = e 2ij(π+Φ) F (γj,j) jζ AP (z A ). (A.44)
The map ζ → Jζ from one family to the other can also be induced using the infinite-dimensional complex structure (A. [START_REF] Martin-Dussaud | A Primer of Group Theory for Loop Quantum Gravity and Spin-foams[END_REF]). An explicit calculation gives in fact

J F (γj,j) jζ (z A ) = e -i arctan γ d j π ||z|| 2(iγj-1-j) [ζ|z 2j = e -i arctan γ F (γj,j) jJζ (z A ), (A.45) J F (γj,j) jJζ (z A ) = (-1) 2j e i arctan γ F (γj,j) jζ (z A ), (A.46)
where we fixed γ > 0 for convenience. When defining the coherent vertex amplitude, it is convenient to have an antipodal map on the coherent state at the target of each link, so that the vectors can be interpreted consistently as outgoing (or incoming, if one prefers) normals. From the previous formulas, we know three different ways to apply the antipodal map: via ζ AP , via J |ζ , or via J F . The resulting amplitudes differ by an overall phase. To keep the notation as simple as possible in the main text, we used the first option, so as to be able to refer only to unit vectors and bypass spinors. The coherent basis elements used in Eq. ( 1) are defined as follows,

D (γj ab ,j ab ) j ab ,-n ab ,j ab , n ba (h -1 a h b ) := D (γj ab ,j ab ) j ab ,j ab ,j ab ,j ab n † (ζ AP ab )h -1 a h b n(ζ ba ) = e -2ij ab Φ ab D (γj ab ,j ab ) j ab ζ ab j ab Jζ ba (h -1 a h b ) = e -2ij ab Φ ab CP 1 dµ(z A ab )F (γj ab ,j ab ) j ab ζ ab (h T a z A ab )F (γj ab ,j ab ) j ab Jζ ba (h T b z A ab ) = e -2ij ab Φ ab d j ab π CP 1 dµ(z A ab ) exp S ab ||h T a z ab || 2 ||h T b z ab || 2 , (A.47)
where in the last step we defined the action At this point we can follow the saddle point approximation developed in [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF], but not before explaining the different conventions used.

S ab := j ab ln h T a z ab | ζab 2 [ ζba |h T b z ab 2 ||h T a z ab || 2 ||h T b z ab || 2 + iγj ab log ||h T b z ab || 2 ||h T a z ab || 2 .
Phase ambiguity and the conventions of [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF] The basis of homogeneous functions used in [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF] differs from Eq. (A.33) in choosing k instead of -k and a different spinorial parametrization,

F B (ρ,k) jm (z A ) := d j π 1 ||z|| 2(1-iρ) D (j) m,k g B (z A ) , g B (z A ) := 1 ||z|| z 0 -z 1 z 1 z0 . (A.51)
In spite of this generic difference, the two choices coincide for the minimal weights used in the EPRL model, since

D (j) m,j g B (z A ) ≡ D (j) m,-j g(z A ) . (A.52)
So there is no difference at this level. The choices however reflect that [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF] takes as fundamental family of coherent states the highest weights, and furthermore the spinorial complex structure (A. [START_REF] Barrett | Quantum gravity asymptotics from the SU(2) 15j symbol[END_REF]) is defined with an opposite sign. The spinorial bilinear in [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF] consistent with (A.45), since we have the same sign convention for J . The link "propagator" is defined in [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF] from the infinite-dimensional bilinear (A.37) (there denoted β(F, F ) but otherwise identically defined). This is the same as a coherent unitary representation matrix up to a phase, which coincides with [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF] (with our definition of J of opposite sign, but this is irrelevant because of the squares). Alternatively, we can reproduce the action of [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF] keeping the same spinors ζ ab = ξ ab and the same group elements h a if we take a global complex conjugate and map γ → -γ.

B EPRL critical point equations and Hessian

We can at this point recall briefly the key steps of the saddle point approximation of [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF], and complete it with the explicit calculation of the Hessian. We look for critical points defined by the vanishing of the gradient of Eq. (A.57) with respect to the variables integrated over, the spinors z ab , and the group elements h a . Vanishing of the spinor gradient gives These are not equations for the integration variables, but directly a restriction on the boundary data. On the other hand, also Eq. (B.59) cannot be solved in general for z ab and h a , but only if the boundary data satisfy special conditions, which leads us to the classification recalled in Sec.2.

In the first subset satisfying Eq. ( 3), let us further specialize the boundary data to the twisted spike configuration [START_REF] Rovelli | Covariant Loop Quantum Gravity[END_REF], where all normals are pairwise antiparallel. Plugging this condition into Eq. (B.59) and looking at the norms, we find equations for the 3D representation matrices H a := D (1) These have a trivial solution with h a = ±1 ∀a. This corresponds in general to a vector geometry in the twisted spike gauge. However if additionally the shape-matching conditions are satisfied, it is possible to find a second, nontrivial solution, given by h (c) 1 = 1 and

h (c) a = ± exp (θ 1a n 1a • σ) (B. 62
)
for a = 1. This is the second critical point for Euclidean Regge data in the twisted spike. As for the spike configuration, one rotates back the tetrahedra so that all edges are aligned, and in that case the two critical points are Eq. ( 20) in the main text. See also Fig. 8 of [START_REF] Donà | SU(2) graph invariants, Regge actions and polytopes[END_REF]. Then to evaluate the action at the critical point, we plug the solutions back in Eq. (B.59) and determine the υ ab phases.

In the second subset, the solutions are not in the SU(2) subgroup, and their derivation is more involved, see [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF]. Again we can simplify the analysis choosing judiciously the boundary data. Fixing the spike configuration described in the main text, and using the freedom to gauge fix h 1 = 1, we can solve the system with Mathematica, finding the set of four solutions [START_REF] Dittrich | Area-angle variables for general relativity[END_REF]. Notice the presence of the iπ: this is due to the inversion required by the fact that the first tetrahedron is past pointing, whereas all the others are future pointing. Plugging them back into Eq. (B.59) we determine the υ ab phases entering the action at the critical point.

For vector geometries in the twisted spike, υ ab = 0. For the Regge geometries with two critical points, the υ ab contribute not only to the Regge action oscillation but also to the overall phase of the amplitude. Their value for the configurations used in this paper are computed in the Mathematica files in [START_REF]The code sl2cfoam is public and available at the[END_REF]. However, the overall phase is the only part of the analytic calculation that does not fit our numerical data. This could be due to an overlooked phase in the previous appendix, or to a mismatch between the analytic definition and the numerical implementation that we were not able to identify. This is of limited harm, since the global phase is irrelevant to the model, and we determined it from a numerical fit for the sake of plotting real data points.

B.1 Hessian

The Hessian H is the matrix of second derivatives of the action (A.57) with respect to the group elements h a and the spinors z ab . It is a 44 × 44 matrix of the form

H = H hh H hz H zh H zz . (B.63)
The generic structure of the Hessian is presented in [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF], however not in a form explicit enough for numerical evaluations. Here we close this gap, and provide the explicit expressions for all the blocks. The top left block is a 24 × 24 block diagonal matrix

H hah b = δ ab H ha ,
where each block is a 6 × 6 matrix given by The off-diagonal block of Eq. (B.63) is a 24 × 20 matrix defined with sparse blocks

H ha = 1 2 b =a j ab - (
H hz =     H h1z12 H h1z13 H h1z14 H h1z15 0 0 0 0 0 0 H h2z12 0 0 0 H h2z23 H h2z24 H h2z25 0 0 0 0 H h3z13 0 0 H h3z23 0 0 H h3z34 H h3z35 0 0 0 H h4z14 0 0 H h4z24 0 H h4z34 0 H h4z45    
where each block is the following 6 × 2 matrix, H hcz ab =δ cb k ab ie i2υ ba v ba 0 iγe i2υ ba v ba (iγ -1) e -i2υ ab v ba + δ ca k ab 0 ie -i2υ ab w ab (iγ + 1) e i2υ ab w ab iγe -i2υ ab w ab defined in terms of the 3D complex vectors By symmetry, the block H zh is a 20 × 24 matrix that can be obtained by transposition of H hz . The expression reported here must be evaluated at the critical point, and since we know that h a h † a is a boost along n 1a , everything can be written as a function of dihedral angles, 3D normals, and phases υ ab . The result is rather cumbersome and we do not write it here, but we refer to the Mathematica notebook SL2Cderivatives in [START_REF]The code sl2cfoam is public and available at the[END_REF] that generates the symbolic expression of the Hessian starting from the action, where this final step is made.

Finally, let us review the various numerical factors that contribute to N c in the asymptotic formula. This includes the inverse square root of the Hessian determinant, as well as the integration measure at the critical point. For the spinor measure Ω ab := ||h † a z ab || -2 ||h † b z ab || -2 dµ(z A ab ), one needs to make a choice of section of the tautological spinor bundle, the same made in evaluating the (inverse square root of the) Hessian at the critical point. Only their product is independent from this choice. The section we chose is ||z ab || = 1 and υ ab = 0 for a < b, and the result can be read from SL2Cderivatives. For the group integrals, we use the normalization of the SL(2, C) Haar measure given in [START_REF] Ruhl | The Lorentz Group and Harmonic Analysis[END_REF]. This gives a factor 1/(2 8 π 4 ) 4 in both the numerical code and the value of the group measure at the critical point. We then have a factor (2π) 22 from the Gaussian integrations, a factor 2 4 from the double multiplicity of the solutions h a = ±h (c) for a = 1, and a factor (2/π) 10 from the normalization of the coherent states in the amplitude. This gives the factor 2 36 π 12 of [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF] times the factor 2 -32 π -16 from the Haar measure. Putting everything together,

N c = 2 4 π 4 1 -det H |σ a<b Ω ab|c . (B.64)
If we revert to our action (A.48), we can derive its Hessian from the one computed in this section by taking the complex conjugate and inverting the sign of γ. Numerical evaluations shows that the determinat of the hessian is invariant under this process. Therefore, the coefficient (B.64) is the same for both (A.48) and (A.57) used in [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF].

Figure 1 :

 1 Figure 1: Classification of the boundary data. It is interesting to compare the situation with SU(2) BF theory [16, 10] and with the Euclidean EPRL model [29, 2, 12]. All three models have the same boundary states and thus some boundary data (j ab , n ab ) in the coherent amplitudes. What changes are the group elements being integrated over [which belong to SL(2, C), SU(2) and SO(4) respectively] and therefore the existence of critical points for a given set. For instance, Lorentzian Regge data only admit critical points for the SL(2, C) model.

{15j} 5 a=2d

 5 ka B γ 4 (j ab , l ab ; i a , k a )

Figure 2 :

 2 Figure 2: Left panel: Parametrization of the tetrahedron with given areas in terms of a dihedral angle φe and the corresponding twist angle τe. Right panel: The dihedral angle φe = φ (a) bc and the corresponding twist angle τe in the Kapovich-Millson polygon made with the scaled normals j ab n ab and jac nac. The variable µe is the length of the dashed line.

Figure 3 :

 3 Figure 3: Left panel: 1 + 1 spacelike Lorentzian simplex and its antispike configuration, obtained boosting two sides in the 1D frame of the base. Two normals are future pointing and one is past pointing. The two hyperbolae are the points you obtain boosting the two sides. Right panel: The same simplex and the spike configuration, obtained acting with a T P transformation after the antispike boosts. All normals are now pointing in the same direction.

Figure 5 :

 5 Figure 5: Convergence of the internal sums at λ = 1 for Euclidean (left panel) and Lorentzian (right panel) data.

Figure 8 :

 8 Figure 8: Study of the impact of truncating three of the five coherent states distributions, for the simplified model with Euclidean equilateral boundary data, and γ = 1.2. Left panel, The error introduced: We compute the relative error σ between the exact amplitude and various truncations corresponding to different values of α -defined in the previous figure, the smaller α the bigger the truncationas a function of the rescaling parameter λ. Setting ourselves an error tolerance of 3%, we truncate at α = 3.Right panel, the improvement on the evaluation times: We compute the ratio between the number of terms Nα=3 surviving in a truncation with α = 3 and the total number Ntot. The impact of the truncation grows roughly linearly in λ (the solid line is 13.6λ -3/2 ), and we see that for instance at λ = 20 only 5% of the terms survive.

Figure 9 :

 9 Figure 9: The truncation relative error in the ∆ = 1 contribution of the shelled summation and Euclidean boundary data.

Figure 10 :

 10 Figure 10: The improvement of the evaluation times obtained truncating one of the Gaussians. Left panel: We perform the same test as shown in Fig.8for the Lorentzian configuration, truncating only on one coherent state (the equilateral one). We can cut only a single tetrahedron due to the fact that, for low rescaling parameters, the isosceles tetrahedra do not approximate a Gaussian, as shown in Fig.7. In this case, the fit gives 0.8λ -1/2 . Right panel: The truncation relative error in the simplified model contribution for Lorentzian boundary data, used in Section 5.4.

Figure 11 :

 11 Figure 11: The amplitude for generic boundary data not admitting critical points, showing an exponential falloff.Both log-log plots showing the absolute value of the amplitude for the simplified model (∆ = 0) and the first shelled sum (∆ = 1), with γ = 1.2. Each plot took roughly 2h for the simplified model, and 36h for ∆ = 1. Left panel: Open data, with spins equal and random normals not satisfying closure. The line is 2.1 × 10 -7 λ -12 , added to help the eye identifying the exponential falloff. Right panel: Closed twisted geometries data, with the normals satisfying closure but not the vector geometry conditions. The line is 3.2 × 10 -12 λ -12 , added to help the eye identifying the exponential falloff.

Figure 12 :

 12 Figure 12: Asymptotics of vector geometries, γ = 1.2.

γN

  c = N P c λ 12 A v 0.1 1.76 × 10 -13 + i1.87 × 10 -14

Figure 17 :

 17 Figure 17: Left panel: Comparison of the simplified model with the EPRL asymptotics and the modified asymptotics

Figure 18 :

 18 Figure 18: Top panel: A pictorial representation of the spike of the Lorentzian 4-simplex used in our numerical analysis. One tetrahedron is equilateral with areas 5, the other four are isosceles with base area 5 and the others 2. Middle panel: Twisted spike for the Euclidean 4-simplices with areas compatible with the spins in the first shell. The 4-simplex on the left has one equilateral tetrahedron with areas 5, one isosceles tetrahedron with base area 5 and the other 3 and three tetrahedra with one area 5, one area 3 and the other two areas 2. The one on the right has one equilateral tetrahedron with areas 5, and four isosceles tetrahedra with base area 5 and the others 3 . Bottom panel: Twisted spike for the Euclidean 4-simplices with areas compatible with the spins in the second shell. The 4-simplex on the left has one equilateral tetrahedron with areas 5, one isosceles tetrahedron with base area 5 and the other 4 and three tetrahedra with one area 5, one area 4 and the other two areas 2. The one in the middle is similar but with areas 5, 4, and 3. The one on the right has one equilateral tetrahedron with areas 5, and four isosceles tetrahedra with base area 5 and the others 4.

29

 29

  ζ| J|j, ζ] = j n, (A.[START_REF]The code sl2cfoam is public and available at the[END_REF] namely the state points in the antipodal direction on the sphere. The two families are in fact related by a parity transformation realized by the complex structureJ (j) |j, ζ := (j) |j, ζ = |j, ζ], J (j)2 = (-1) 2j . mn := (-1) j-m δ m,-n = D (j)mn (iσ 2 ) is the spinorial metric. Another useful operation is the antipodal map

(

  γj ab ,j ab ) j ab ,m ab ,j ab ,m ba (h -1 a h b ). (A.41)

(A. 48 )(j ab π 5 a=2 dh a CP 1

 4851 Therefore our definition (1) of the coherent vertex amplitude isA v (j ab , n ab ) = γj ab ,j ab ) j ab ,-n ab ,j ab , n ba (h -1 a h b ) = = e -2i a<b j ab Φ ab a<b d a<b dµ(z A ab ) ||h T a z ab || 2 ||h T b z ab || 2 exp S,

J

  is [w, z] = [z|w .The coherent basis is defined from the highest weights, and denotedIφ ξ (z A ) := d j π ||z|| 2(iγj-1-j) [ξ|z 2j ≡ F Iφ ξ (z A ) = (-1) 2j e -i arctan γ d j π||z|| 2(iγj-1-j) ξ|z 2j ≡ (-1) 2j e -i arctan γ F (γj,j) jξ (z A ), (A.54)

P 2 + 2 ,

 22 ab := (h a Iφ ab , h b Iφ ba ) ≡ (h a J Iφ ab , h b Iφ ba ) = (-1) 2j ab e i arctan γ D (γj ab ,j ab ) j ab ξab j ab Jξ ba(h -1 a h b ) = (-1) 2j ab e i arctan γ CP 1 dµ(z A )F (γj ab ,j ab ) j ab ξab (h T a z A ) F (γj ab ,j ab ) j ab Jξ ba (h T b z A ). (A.55)Hence our expression (A.47) for the link propagator differs from the one of[START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF] by the overall phase exp{i a<b (2j ab (π + Φ ab ) + arctan γ)} reported in footnote 4, and by the identificationζ ab = ξab . (A.56) With this map, and the renaming h a → ha , the action becomes S ab := j ab ln h † a z ab |ξ ab 2 Jξ ba |h † b z ab 2 ||h † a z ab || 2 ||h † b z ab || iγj ab log ||h † b z ab || 2 ||h † a z ab ||

|ξ ab = e iυ ab ||h † a z ab || h † a |z ab , |Jξ ba = e iυ ba ||h † b z ab || h † b |z ab , (B. 58 )

 58 which can be combined to(h † a ) -1 |ξ ab = ||h † b z ab || ||h † a z ab || e i(υ ab -υ ba ) (h † b ) -1 |Jξ ba , (B.59a)as well ash a |ξ ab = ||h † a z ab || ||h † b z ab || e i(υ ab -υ ba ) h b |Jξ ba . (B.59b)Here υ ab are phases to be determined. Because of the homogeneity of the action, these are only four independent real equations per link. Two fix z ab up to the irrelevant choice of section, and the remaining ones fix the group element and υ ab phases. On-shell of these equations, the vanishing of the h a gradient gives the closure conditions b =a j ab n ab = 0. (B.60)

  (h a ),H a n ab = -H b n ba = H b n ab . (B.61)

2 +|ξ ba 2 ξ

 22 j ab (iγ -1) e i2(2υ ba +arctan γ) Jξ ba |h b h † b ba |(h b h † b ) 2 |ξ ab + | ξ ba |h b h † b |ξ ab | 2 ξ ba |(h b h † b ) 2 |ξ ab + | ξ ba |h b h † b |ξ ab | 2 -e -i4υ ba ξ ba |h b h † b |Jξ ba 2 .

v

  ba = Jξ ba | σh b h † b |ξ ba -n ba Jξ ba |h b h † b |ξba , and w ab = Jξ ab |h a h † a σ|ξ ab + n ab Jξ ab |h a h † a |ξ ab .

  : One gives the ten areas, one dihedral angle φ

	between two normals belonging to different tetrahedra; the latter angle is not gauge invariant and has to
	satisfy an inequality. It is convenient to choose the gauge (4) with all normals pairwise antiparallel. A
	simple configuration following this construction, with all areas equal and not giving a Regge geometry, is
	the following:
		φ	(2) 34 =	6 15	π,	φ (3) 25 =	8 15	π,	φ (4) 25 =	7 15	π,	φ (5) 34 =	8 15	π,	φ 53,42 =	9 15	π,	(11)
	where φ 53,42 := arccos ( n 53 • n 42 ) is the gauge-dependent angle. The resulting normals are summarized in
	Table 3.
	ñ b b		1				2				3				4			5
	1					(0.82, -0.45, -0.35) (-0.96, 0.05, 0.28)	(-0.76, 0.54, -0.34)	(0.9, -0.15, 0.41)
	2	(-0.82, 0.45, 0.35)				(0.31, 0.21, -0.93)			(1, 0, 0)		(-0.49, -0.66, 0.57)
	3	(0.96, -0.05, -0.28) (-0.31, -0.21, 0.93)					(-0.34, -0.69, -0.64)	(-0.31, 0.95, 0)
	4	(0.76, -0.54, 0.34)	(-1, 0, 0)	(0.34, 0.69, 0.64)					(-0.1, -0.14, -0.98)
	5	(-0.9, 0.15, -0.41)	(0.49, 0.66, -0.57)	(0.31, -0.95, 0)			(0.1, 0.14, 0.98)
									(a) bc in four different tetrahedra, plus one additional angle

Table 4 :

 4 Choice of normals for the Euclidean Regge data, corresponding to the twisted spike. These are reconstructed from an equilateral 4-simplex using the algorithm described in the text.

  below.

	ñ b b	1	2	3	4	5
	1		(1, 0, 0)	(-0.33, 0.94, 0)	(-0.33, -0.47, 0.82)	(-0.33, -0.47, -0.82)
	2	(-1, 0, 0)		(0.83, 0.55, 0)	(0.83, -0.28, 0.48)	(0.83, -0.28, -0.48)
	3	(0.33, -0.94, 0)	(0.24, 0.97, 0)		(-0.54, 0.69, 0.48)	(-0.54, 0.69, -0.48)
	4	(0.33, 0.47, -0.82)	(0.24, -0.48, 0.84)	(-0.54, 0.068, 0.84)		(-0.54, -0.76, 0.36)
	5	(0.33, 0.47, 0.82)	(0.24, -0.48, -0.84) (-0.54, 0.068, -0.84) (-0.54, -0.76, -0.36)

Table 6 :

 6 3.53 × 10 -13 cos 0.963 × 10 -2 λ + 0.106 1.2 2.82 × 10 -16 + i5.7 × 10 -16 1.27 × 10 -15 cos (0.116λ + 1.111) 7 8.86 × 10 -25 -i3.67 × 10 -25 1.92 × 10 -24 cos (0.674λ -0.392) Asymptotics for Lorentzian boundary data at different values of γ.

	8	1.48 × 10 -25 -i4.32 × 10 -26	3.08 × 10 -25 cos (0.770λ -0.284)

  1 -n ab ⊗ n ab ) -i ab (1 -n ab ⊗ n ab ) + ( ab + iγ) n ab -i ab (1 -n ab ⊗ n ab ) -( ab + iγ) n ab -(1 + 2i ab γ) (1 -n ab ⊗ n ab )where ab = sign (b -a) and is the Hodge dual in the 3D internal space. The bottom right block of Eq. (B.63) is a 20 × 20 block diagonal matrix H z ab z cd = δ ac δ bd H z ab z ab where each block is a 2 × 2 matrix H z ab =j ab (iγ + 1) e i4υ ab ξ ab |h a h † a |Jξ ab 2 | Jξ ab |h a h † a |Jξ ab | 2 -Jξ ab |(h a h † a ) 2 |Jξ ab | Jξ ab |h a h † a |Jξ ab | 2 -Jξ ab |(h a h † a ) 2 |Jξ ab -e -2i(2υ ba +arctan γ) Jξ ab |h a h † a |ξ ab

To be more precise, twisted geometries have additional 'twist' angle variables, denoted ξ ab and canonically conjugated to the spins[21,

[START_REF] Freidel | From twistors to twisted geometries[END_REF]]. These angles do not show up in the coherent amplitude[START_REF] Engle | LQG vertex with finite Immirzi parameter[END_REF], defined by eigenvectors of the Casimir spin operator. They would appear in amplitudes where also the areas are semiclassical and not quantum, like in the graviton propagator calculations, see e.g.[START_REF] Bianchi | Graviton propagator in loop quantum gravity[END_REF][START_REF] Speziale | Background-free propagation in loop quantum gravity[END_REF]. The boundary data of (1) are thus a subset of twisted geometries with ξ ab = 0.2 That this is compatible with Eq. (3) is not obvious, and one of the key results of the reconstruction.

As for full Lorentz invariance, this is broken by the Y -map identifying k = j. Nonetheless, it can be shown that the vertex amplitude transforms covariantly under Lorentz transformations of the boundary data, and when various vertex amplitudes are glued together to form an extended spin foam, Lorentz invariance is restored in the bulk[START_REF] Rovelli | Lorentz covariance of loop quantum gravity[END_REF].

The definition of the coherent vertex amplitude used here differs from the one in[START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF] by an overall phase, given for γ > 0 by exp{i a<b (2j ab (π + Φ ab ) + arctan γ)}. This is due to a phase difference between the Wigner matrices and the antilinear pairing of two unitary representations, as well as to our use of the antipodal map instead of the complex structure in associating coherent states to the target nodes. See Appendix A for details.

As was done for the graviton propagator with the Euclidean Barrett-Crane model[START_REF] Christensen | Numerical evidence of regularized correlations in spin foam gravity[END_REF].

The expression of[START_REF] Collet | A (simple) expression of the unitary-irreducible SL(2,C) representations as a finite sum of exponentials[END_REF] is very similar to the one previously derived in[START_REF] Rashid | Boost Matrix Elements Of The Homogeneous Lorentz Group[END_REF], but with the crucial difference that it only has finite sums. It is obtained starting from the integral representation, whereas the one of[START_REF] Rashid | Boost Matrix Elements Of The Homogeneous Lorentz Group[END_REF] starts from the infinite power series for the hypergeometric functions.

Note that the resulting function from vertices (or equivalently from lengths) to areas is always injective and generically bijective; it fails to be bijective in special cases, like precisely for all equal areas, that can correspond to either an equilateral 4-simplex, or a 4-simplex with seven equilateral triangles and three isosceles. The 3D normals then distinguish the two cases.

This is determined computing the Euclidean scalar product between the 4D normal and the vector connecting the chosen vertex to the center of the 4-simplex. Throughout the paper, we always work with outgoing normals, both 3D and 4D, and accordingly all dihedral angles are external.

A tetrahedron on the boundary of a Lorentzian 4-simplex is spacelike if its four-dimensional normal Na is timelike, N I a N aI = -1. Conversely, the tetrahedron is said to be timelike if its four normal is spacelike.

More precisely slightly after, at ī = (2/ √ 3)j 1.15j[START_REF] Livine | A New spinfoam vertex for quantum gravity[END_REF].

The reader may wonder why we can use open configurations in the boundary states, since these represent coherent intertwiners and closure is the classical counterpart of the quantum gauge invariance. The answer is that the normals are just classical labels, therefore one can have gauge invariant states with nonclosed normals. These states are on the other hand over-redundant, and can be effectively removed without losing coherence. The result is a family of coherent states depending only on the cross ratio parametrizing the intrinsic shape of the tetrahedron[START_REF] Freidel | Holomorphic Factorization for a Quantum Tetrahedron[END_REF]. This family was used for numerical calculations of the volume operator in[START_REF] Bianchi | Polyhedra in loop quantum gravity[END_REF], but we refrain from using it here to keep the treatment as simple as possible.

The two critical points are related by a parity transformation, hence the notation Pc. Then one finds Sc = iλ(S R + Φc) and SPc = iλ(-S R + Φc), leading to Eq.[START_REF] Rashid | Boost Matrix Elements Of The Homogeneous Lorentz Group[END_REF].

More precisely, one gets angle-matching conditions. However, for triangular faces with the areas already matching by construction, matching angles imply matching shapes. This is not the case for general polytopes with nontriangular faces, and in fact the large spin asymptotics of nonsimplicial vertices contains not just Regge geometries, but more general conformal twisted geometries, see the discussion in[START_REF] Donà | SU(2) graph invariants, Regge actions and polytopes[END_REF].

A word of caution to avoid possible confusion: If one plots the asymptotic formula as a real function, the real frequency given by S R is much higher than the frequency of oscillations that can be deduced interpolating the (half-)integer sampling, a situation familiar from the study of SU(2) asymptotics. Hence when we speak of matching frequency of oscillations between the numerical data and the analytic formula, we refer to the "effective" frequency of the (half-)integer oscillations, which is not the Regge cation, but a function thereof.

This is achieved through a further imposition of the Y map inside the product h -1 a h b in Eq. (6), see[START_REF] Speziale | Boosting Wigner's nj-symbols[END_REF] for details and motivations.

With a peak when the virtual ρ label takes value γj, thus showing a certain type of simplicity off-shellness in the model.

The complex conjugation in z|, which could be avoided by defining the coherent states as m F (ρ,j) jm (z A ) j, ζ|j, m , is actually all right, and familiar from (holomorphic) representations generated by (unnormalized) coherent states, e.g. fn(z) := z|n = z n /n! for the harmonic oscillator.

Acknowledgments

We would like to thank Jean-Roch Liebgott and Lorenzo Bosi for help with the server at CPT, and Francesco Gozzini and François Collet for discussions and the shared work on the numerical code. P.D. thanks Tommaso De Lorenzo and Eugenio Bianchi for useful discussions, and support from FQXi grants No. 2016-165616 and No. 2018-190485. P.D. is also supported by the NSF grants No. PHY-1505411 and No. PHY-1806356, and the Eberly research funds of Pennsylvania State University. We thank Muxin Han, Hongguang Liu and Dongxue Qu for exchanges on the calculation of the Hessian that led us to correct the sign of its phase.

A Conventions and explicit formulae

In this appendix we fix our conventions and notation, and provide the explicit relation between our version [START_REF] Dona | Numerical methods for EPRL spin foam transition amplitudes and Lorentzian recoupling theory[END_REF] of the coherent amplitude and the one in [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF]. The difference is just in the phase and a complex conjugation of the spinors, but it has multiple origins due to the many conventional choices involved. Our conventions follow those of [START_REF] Perelomov | Generalized coherent states and their applications[END_REF][START_REF] Ruhl | The Lorentz Group and Harmonic Analysis[END_REF][START_REF] Varshalovich | Quantum Theory of Angular Momentum: Irreducible Tensors, Spherical Harmonics, Vector Coupling Coefficients, 3nj Symbols[END_REF].

A.1 SU(2) coherent states and spinors

SU(2) coherent states

The uncertainty in the direction of a state |j, m in the SU(2) irrep j is minimized by the lowest and highest weights m = ±j. Acting with a rotation g ∈ SU(2) on |j, ±j defines two alternative families of coherent states. Orbits of the isotropy subgroup U(1) ⊂ SU(2) generate an irrelevant global phase in the coherent states, and it is sufficient to consider two-parameters group elements in the form [START_REF] Perelomov | Generalized coherent states and their applications[END_REF]

This is also known as the Hopf section for the SU(2) S 2 × S 1 fibration, and we parametrized ζ as the stereographical projection of the sphere from the south pole, with (Θ, Φ) the zenith and azimuth angles. In terms of Wigner's matrices as D(α, β, γ) = e -iαJz e -iβJy e -iγJz [START_REF] Varshalovich | Quantum Theory of Angular Momentum: Irreducible Tensors, Spherical Harmonics, Vector Coupling Coefficients, 3nj Symbols[END_REF] (caveat in reading the files in [START_REF]The code sl2cfoam is public and available at the[END_REF], Wolfram's Mathematica uses opposite sign conventions for this parametrization),

The standard choice is to define coherent states starting from the lowest weight [START_REF] Perelomov | Generalized coherent states and their applications[END_REF],

They satisfy j, ζ| J|j, ζ = -j n, n := (sin Θ cos Φ, sin Θ sin Φ, cos Θ), (A.4)

and the important factorization property

The phase ambiguity associated with the U(1) subgroup, namely with choosing a different section than Hopf's, is the one that can be used to make the coherent amplitude real for any boundary data configuration, as discussed in [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF] and in the main text.

It is often convenient to use also the second family, the one constructed from the highest weight. To distinguish it, we use a different ket notation,

Spinors

The fundamental irrep j = 1/2 is isomorphic to C 2 , and the relation of the coherent states to spinors is very useful to study properties of the SL(2, C) unitary irreps. We denote a spinor with z A ∈ C 2 , A = 0, 1, and use

to raise and lower indices. We equip C 2 with a complex structure J and define the dual spinor

Notice that so defined, the complex structure is consistent with Eq. (A.8) for j = 1/2. AB gives an SL(2, C)-invariant bilinear coupling, and we can introduce an Hermitian norm given by complex conjugation (book kept by dotting the indices) and the identity matrix δ A Ȧ, which is nothing but a choice of canonical timelike vector in Penrose's abstract index notation. With a Dirac-like notation for spinors [START_REF] Speziale | The twistorial structure of loop-gravity transition amplitudes[END_REF], the relations between spinors and their dual can be compactly written

with = iσ 2 , and the Hermitian and Lorentz-invariant bilinear are

Finally, for the norm we use

Coming back to the SU(2) coherent states, the lowest-weight ones define spinors

of unit norm and arg z 1 = 0. The map z 0 = -sin Θ 2 e -iΦ , z 1 = cos Θ 2 provides a choice of section for the C 2 CP 1 × CP 1 fibration. Conversely, an arbitrary spinor defines a coherent state for the fundamental irrep which is redundantly parametrized by an additional norm and phase,

The additional phase corresponds to the phase ambiguity in the coherent states previously discussed, whereas the norm can be eliminated by working always with normalized spinors. Spinorial SU(2) coherent states feature prominently in spin foam models and their applications, see e.g. [START_REF] Dupuis | Holomorphic Simplicity Constraints for 4D Spinfoam Models[END_REF][START_REF] Dupuis | Holomorphic Lorentzian Simplicity Constraints[END_REF][START_REF] Speziale | The twistorial structure of loop-gravity transition amplitudes[END_REF][START_REF] Hnybida | Spin Foams Without Spins[END_REF][START_REF] Banburski | Pachner moves in a 4D Riemannian holomorphic Spin Foam model[END_REF][START_REF] Bonzom | Duality between Spin networks and the 2D Ising model[END_REF]]. An immediate advantage of their use is to write the resolution of the identity as Gaussian integrals.

For the highest weight,

Again, an arbitrary dual spinor defines a redundant coherent state with an additional norm and phase,

The expressions (A.21) and (A.23) can be generalized to an arbitrary irrep j, to provide spinorial coherent states for SU [START_REF] Freidel | A New Spin Foam Model for 4D Gravity[END_REF]. To that end, we consider the matrix

1). (A.24)

We then define the spinorial coherent states

where in the second equality we used the following properties of the Wigner matrices:

(A.27)

These expressions reduce to Eqs. (A.21) and (A.23) for j = 1/2. Finally, we note the factorization properties

A.2 SL(2, C) unitary irreps of the principal series

Homogeneous representation

We follow the conventions of [START_REF] Ruhl | The Lorentz Group and Harmonic Analysis[END_REF], with the notational difference

Unitary irreps of the principal series of SL(2, C) are labeled by a pair (ρ ∈ R, k ∈ Z/2) [START_REF] Gelfand | Representations of the rotation and Lorentz groups and their applications[END_REF][START_REF] Naimark | Linear representations of the Lorentz group[END_REF][START_REF] Ruhl | The Lorentz Group and Harmonic Analysis[END_REF], and can be represented on a space of homogeneous functions of two complex variables,

The scalar product is defined choosing any section of the

The scaling of the Lorentz-invariant measure dµ(z A ) guarantees that the whole integrand is homogeneous of degree zero, and thus the integration is independent of the choice of section.

For the EPRL model with spacelike tetrahedra we pick Naimark's orthonormal basis, labeled by the eigenfunctions of L 2 and L z and given in the conventions of [START_REF] Ruhl | The Lorentz Group and Harmonic Analysis[END_REF] by

Here D (j) are Wigner matrices, and g(z A ) is the same matrix defined in Eq. (A.24). The function Ψ ρ j parametrizes a freedom in choosing the phase of the basis elements. This is set to zero in [START_REF] Ruhl | The Lorentz Group and Harmonic Analysis[END_REF] and [START_REF] Barrett | Lorentzian spin foam amplitudes: Graphical calculus and asymptotics[END_REF]. This economical choice however leads to complex Clebsch-Gordan coefficients. An alternative choice is to take [START_REF] Speziale | Boosting Wigner's nj-symbols[END_REF] e iψ ρ j = (-1) -j 2 Γ(j + iρ + 1) |Γ(j + iρ + 1)| , (A. [START_REF] Rashid | Boost Matrix Elements Of The Homogeneous Lorentz Group[END_REF] which guarantees reality of the Clebsch-Gordan coefficients, and as a consequence also of all EPRL amplitudes A v (j ab , i a ) in the orthogonal intertwiner basis. This is the convention used in our numerical code, but this phase difference goes away in the expression of the coherent amplitude.