Pietro Donà

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

Numerical methods for EPRL spin foam transition amplitudes and Lorentzian recoupling theory

The intricated combinatorial structure and the non-compactness of the Lorentz group have always made the computation of SL(2, C) EPRL spin foam transition amplitudes a very hard and resource demanding task. With sl2cfoam we provide a C-coded library for the evaluation of the Lorentzian EPRL vertex amplitude. We provide a tool to compute the Lorentzian EPRL 4-simplex vertex amplitude in the intertwiner basis and some utilities to evaluate SU(2) invariants, booster functions and SL(2, C) Clebsch-Gordan coefficients. We discuss the data storage, parallelizations, time, and memory performances and possible future developments.

Introduction

The spin foam framework provides an approach to the dynamics of loop quantum gravity that is both background independent and Lorentz covariant. The most successful model is the one proposed by Engle, Pereira, Rovelli and Livine (EPRL) [START_REF] Engle | LQG vertex with finite Immirzi parameter[END_REF][START_REF] Etera | A New spinfoam vertex for quantum gravity[END_REF][START_REF] Etera | Consistently Solving the Simplicity Constraints for Spinfoam Quantum Gravity[END_REF] and independently by Freidel and Krasnov (-FK) [START_REF] Freidel | A New Spin Foam Model for 4d Gravity[END_REF]. The model is formulated as a sum over simplicial triangulations with given boundary providing transition amplitudes for (projected) spin network states of the canonical theory restricted to four valent nodes. This procedure is equivalent to a sum over histories of quantum geometries providing in this way a regularised version of the quantum gravity path integral.

A generalization to arbitrary complexes and spin-networks has been developed [START_REF] Kaminski | Spin-Foams for All Loop Quantum Gravity[END_REF]. Nevertheless, in this work, we focus primarily on the simplicial formulation, and we provide all the tools for more general computations.

Computations in this framework are particularly complex. Two of the most successful analytical calculations are: On the one hand, the presence of the discrete Regge action in the spin foam amplitude of a semiclassical state [START_REF] Barrett | Asymptotic analysis of the eprl four-simplex amplitude[END_REF][START_REF] Barrett | Asymptotic analysis of Lorentzian spin foam models[END_REF]. Regge calculus provides a discretization of general relativity based on a triangulation. 1 On the other hand, the graviton n-point function reproduces at the leading order the one of Regge calculus in the semi-classical limit [START_REF] Bianchi | Graviton propagator in loop quantum gravity[END_REF][START_REF] Speziale | Background-free propagation in loop quantum gravity[END_REF][START_REF] Bianchi | LQG propagator from the new spin foams[END_REF][START_REF] Bianchi | Lorentzian spinfoam propagator[END_REF]. This result is of particular historical importance since it revealed the shortcomings of the Barrett-Crane model and led to the development of the EPRL-FK model. In particular, it has been shown that within the Barrett-Crane model it was not possible to reproduce the correct long-distance scaling of the graviton two-point function [START_REF] Alesci | The Complete LQG propagator. I. Difficulties with the Barrett-Crane vertex[END_REF][START_REF] Alesci | The Complete LQG propagator. II. Asymptotic behavior of the vertex[END_REF][START_REF] Alesci | LQG propagator: III. The New vertex[END_REF]. Therefore, the model is unable to implement the simplicity constraints correctly. These are responsible for reducing the topological BF theory, the starting point of any spin foam model, to gravity.

Evaluating spin foam transition amplitudes for assigned boundary data is an arduous task. All the extraordinary results we mentioned are asymptotic computations performed in the large spin limit. At the present time, numerical calculations within the theory are very rare. Spin correlations have been computed with the old Barrett-Crane model [START_REF] Daniel Christensen | Numerical evidence of regularized correlations in spin foam gravity[END_REF][START_REF] Daniel Christensen | Sub-leading asymptotic behaviour of area correlations in the Barrett-Crane model[END_REF] or with the Lorentzian EPRL model but restricted to simple generalized spin foams [START_REF] Sarno | 2-vertex Lorentzian Spin Foam Amplitudes for Dipole Transitions[END_REF].

In this work, we introduce the library sl2cfoam where we provide the tools to evaluate the EPRL spin foam vertex amplitude numerically. It is based on the decomposition of the EPRL vertex amplitude in terms of booster functions introduced in [START_REF] Speziale | Boosting wigner's nj-symbols[END_REF]. These utilities are indicated for the study of vertex amplitudes with coherent boundary data and non, opening new windows on the covariant formulation of loop quantum gravity.

Some preliminary results have already been achieved. Using a primitive version of this code, it was possible to verify, for the first time, the asymptotic behavior of the topological BF SU(2) spin foam amplitude [START_REF] Dona | SU(2) graph invariants, Regge actions and polytopes[END_REF], but also to find the semiclassical region at relatively small values of the spins. The transition amplitudes and spin correlations between two dipole spin networks, in the Lorentzian EPRL model, with up to two (nonsimplicial) vertices has been studied numerically in [START_REF] Sarno | 2-vertex Lorentzian Spin Foam Amplitudes for Dipole Transitions[END_REF]. Even if restricted to simple vertices this was the first full evaluation of an EPRL transition amplitude. The numerical computation of the non-isotropic scaling of the booster functions, performed with a piece of this software, is a central ingredient in the recent estimation of the degree of divergence of the EPRL model [START_REF] Dona | Infrared divergences in the EPRL-FK Spin Foam model[END_REF].

Moreover, we are finishing the numerical exploration of the Euclidean sector of the asymptotic of the Lorentzian EPRL vertex amplitude, reproducing with great accuracy the analytic asymptotic formula [START_REF] Dona | Numerical studies of the Lorentzian EPRL vertex amplitude[END_REF].

With the access to a reliable and systematic way to perform numerical computation within the Lorentzian spin foam framework, in particular with more than a single vertex (or equivalently with many simplices), we are opening a window on a plethora of open questions. One question we plan to address soon is the so-called flatness problem, firstly mentioned by Freidel and Conrady [START_REF] Conrady | On the semiclassical limit of 4d spin foam models[END_REF], later by Bonzom [START_REF] Bonzom | Spin foam models for quantum gravity from lattice path integrals[END_REF] and Hellmann and Kaminski [START_REF] Hellmann | Holonomy spin foam models: Asymptotic geometry of the partition function[END_REF]. They argue that the classical geometries reproduced by the EPRL model in the semiclassical limit are restricted to be flat and, as a result, the model itself is put in question, since GR in four dimensions admits curved spacetime solutions. A preliminary answer, in the Euclidean setting, has been given [START_REF] José | EPRL/FK Asymptotics and the Flatness Problem[END_REF][START_REF] Vincent Bayle | Short-scale Emergence of Classical Geometry[END_REF], showing that curved geometric configurations are asymptotically allowed. We think that a numerical exploration is viable in Lorentzian setting and could help us shedding light on this issue.

The realm of applications of sl2cfoam is not limited to technical results. The computation of physical observables requires the ability to evaluate spin foam transition amplitude with multiple vertices. The study of the black hole to white hole tunneling processes [START_REF] Christodoulou | Planck star tunneling time: An astrophysically relevant observable from background-free quantum gravity[END_REF][START_REF] Christodoulou | Characteristic Time Scales for the Geometry Transition of a Black Hole to a White Hole from Spinfoams[END_REF] and quantum cosmological models [START_REF] Bianchi | Towards Spinfoam Cosmology[END_REF][START_REF] Vidotto | Many-nodes/many-links spinfoam: the homogeneous and isotropic case[END_REF][START_REF] Sarno | 2-vertex Lorentzian Spin Foam Amplitudes for Dipole Transitions[END_REF] is the beginning of a long term project that aims to connect the quantum theory of gravity with observation. We hope that the release of the tools in sl2cfoam will help these computations to flourish.

The purpose of the present work goes beyond giving the tools for the computation of Lorentzian EPRL transition amplitudes. In particular, we believe that the techniques we introduced and the algorithms we implemented can be used by the loop quantum gravity community on problems outside the spin foam context. For example, the numerical evaluation of SL(2, C) Clebsch Gordan coefficients is a challenging and interesting mathematical problem on its own and the very efficient way we introduced to deal with SU(2) invariant can play a key role in exploring numerically the action of the hamiltonian operator in the canonical approach [START_REF] Alesci | Matrix elements of Lorentzian Hamiltonian constraint in loop quantum gravity[END_REF]. We assume that the reader is familiar with the EPRL model, and refer to the original literature [START_REF] Engle | LQG vertex with finite Immirzi parameter[END_REF][START_REF] Etera | A New spinfoam vertex for quantum gravity[END_REF][START_REF] Etera | Consistently Solving the Simplicity Constraints for Spinfoam Quantum Gravity[END_REF][START_REF] Freidel | A New Spin Foam Model for 4d Gravity[END_REF]] and existing reviews (e.g. [START_REF] Perez | The spin foam approach to quantum gravity[END_REF][START_REF] Rovelli | Covariant loop quantum gravity: an elementary introduction to quantum gravity and spinfoam theory[END_REF]) for motivations, details and its relation to loop quantum gravity. We give a short review of the transition amplitude in Section 2.

In Section 3 we present in details how the library sl2cfoam is structured. In Section 3.1 we give a general overview of the problems we faced and the solutions we propose, while in Sections 3.2, 3.3 and 3.4 we discuss its three main components, respectively the evaluation of SU(2) invariants, booster functions and finally the complete vertex amplitude.

In the last two sections we present two additional utilities that computes the Livine-Speziale coherent states in Section 4.1 and the SL(2, C) Clebsch-Gordan coefficients in Section 4.2.

All the tests shown in this paper have been performed on Intel(R) Xeon(R) E5-2687W v2 @ 3.40GHz CPUs provided by CPT servers in Marseille. If not stated otherwise the tests are performed on a single core.

The library sl2cfoam can be found at: https://bitbucket.org/giorgiosarno/sl2cfoam-1.0/. We provide a README file in which we explain how to use the code. We include the test programs we used to perform the benchmarks in this paper and can be seen as examples. The computation of a vertex amplitude for given data has to be computed only once. We include a compressed table of pre-calculated symbols, building blocks, and full amplitudes that we already computed in the repository. We will continuously update it in the future.

The EPRL spin foam transition amplitude

The spin foam formalism introduces a partition function of a given closed 2-complex C as a state sum over SU (2) spins j f on the faces and intertwiners i e on the edges:

Z C = j f ,ie f A f (j f ) e A e (i e ) v A v (j f , i e ) . (1)
defined in terms of the face amplitude A f , and the edge amplitude A e and the vertex amplitude A v . Requiring the correct convolution property of the path integral at fixed boundary, the form of the face amplitude A f (j f ) := 2j f + 1 and the edge amplitude A e (i e ) = 2i e + 1 are fixed [START_REF] Bianchi | Face amplitude of spinfoam quantum gravity[END_REF]. In this paper, we will focus on the numerical evaluation of the vertex amplitude A v defined by the EPRL model. The vertex amplitude is built from the topological SL(2, C) BF spin foam vertex amplitude by imposing, weakly, the simplicity constraints resulting into a restriction of the unitary irreducible representations in the principal series [START_REF] Engle | LQG vertex with finite Immirzi parameter[END_REF][START_REF] Engle | The Loop-quantum-gravity vertex-amplitude[END_REF]. To evaluate it in its original form one should perform four group integrals (one of the original five is redundant and has to be removed to guaratee finiteness [START_REF] Engle | Regularization and finiteness of the Lorentzian LQG vertices[END_REF]). Each group integral is, in general, a six dimensional unbounded highly oscillating integral for which numerical integration methods are very inefficient. To solve this issue, a more computational friendly form for the amplitude, has been proposed [START_REF] Speziale | Boosting wigner's nj-symbols[END_REF]:

A v (j f , i e ) = l f ,ke e (2k e + 1) B 4 (j f , l f ; i e , k e ) {15j} j f ,i4 (l f , k e ) (2) = l f ,ke k 1 k 2 k 5 i 4 k 3 j 2 j 4 j 10 j 9 i 1 i 2 i 5 i 3 j 7 j 1 j 5 j 3 j 2 j 8 j 7 j 6 j 4 j 8 j 3 j 1 j 6 j 5 j 10 j 9
, resulting in a superposition of SU (2) {15j} symbols weighted by one booster functions B 4 per edge e in the considered vertex. To obtain this form each SL(2, C) integral is decomposed in a SU(2) "rotation", a boost along a fiducial axes and another SU(2) "rotation". The compact integrals are then evaluated exactly composing the SU(2) invatiants in [START_REF] Etera | A New spinfoam vertex for quantum gravity[END_REF].

The summation is over a set of auxiliary spins l f2 for each face involving the vertex but excluding the gauge fixed edge for a total of 6 distinct l f , with lower bound l f ≥ j f , and a set of auxiliary intertwiners k e for each edge in the vertex excluding the gauge fixed one for a total of 4 that can assume all the values compatible with triangular inequalities.

The expression (2) is far more convenient than the original one for numerical evaluation. We decomposed the amplitude into smaller building blocks, which are easier to compute. We trade four six dimensional unbounded integrals with a summation of SU(2) invariants, which are widely studied, and a class of onedimensional integrals that we need to perform and we can characterize independently.

These formulas can be generalized to arbirtary valency, for more details see [START_REF] Sarno | 2-vertex Lorentzian Spin Foam Amplitudes for Dipole Transitions[END_REF]. In the library sl2cfoam , in addition to the simplicial amplitude, we also provide the tools to evaluate amplitudes of four valent vertices. The reader can be unfamiliar with the B 4 functions, these are the non compact remants of the SL(2, C) integrals and they encode all the details of the EPRL model. We named them booster functions and are defined in the following way:

B 4 (j f , l f ; i, k) = 1 4π pa j f p f (i) ∞ 0 dr sinh 2 r n a=1 d (γj f ,j f ) j f l f p f (r) l f p f (k) = j 4 i j 3 j 2 j 1 l 4 k l 3 l 2 l 1 , (3)
where d (ρ,k) (r) are the boost matrix elements for γ-simple irreducible representations of SL(2, C) in the principal series, γ is the Immirzi parameter and the (4jm) symbols are reported in Appendix A. All the test and computations in this paper are performed with a conventional value of γ = 1.2. The explicit form of the boost matrix elements can be found in the literature in its general form [START_REF] Ruhl | The Lorentz group and harmonic analysis[END_REF][START_REF] Naimark | Linear Representations of the Lorentz Group[END_REF][START_REF] Rashid | Boost Matrix Elements of the Homogeneous Loretz Group[END_REF][START_REF] Speziale | Boosting wigner's nj-symbols[END_REF], here we just report them in the case of simple irreducible representation:

d (γj,j) jlp (r) =(-1) j-l 2 Γ (j + iγj + 1) |Γ (j + iγj + 1)| Γ (l -iγj + 1) |Γ (l -iγj + 1)| √ 2j + 1 √ 2l + 1 (j + l + 1)! (2j)!(l + j)!(l -j)! (l + p)!(l -p)! (j + p)!(j -p)! 1/2 × e -(j-iγj+p+1)r s (-1) s e -2sr s!(l -j -s)! 2 F 1 [l + 1 -iγj, j + p + 1 + s, j + l + 2, 1 -e -2r ] . (4)
The study of the booster function is still in the initial stage, more work is necessary to fully understand all their properties. For example, while the symmetries are known, the recursion relations still need to be explored. The study of their semiclassical limit is in a preliminary stage [START_REF] Fanizza | Asymptotics of SL(2,C) tensor invariants[END_REF] but an appealing geometrical picture seems to emerge. On one hand, the introduction of booster functions brings a drastic simplification to the computation of the vertex amplitude because it reduces the problem of dealing with many high oscillatory integrals to the study a family of one dimensional integrals, which are easier to handle and manipulate. On the other hand, we introduce a large amount of convergent summation in the computation which need to be dealt with. The net gain of this operation is evident once we realize that each element of the sum can be computed in a fast way and with high precision.

3 The sl2cfoam library

General strategy

We can summarize the evaluation of the vertex amplitude (2) in three main ingredients that we need to mix in the right way:

• all the necessary {15j} symbols,

• all the required booster functions,

• the sum over the auxiliary spins and intertwiners.

In the next sections we will discuss in great detail the problems we encounter in these three steps and their possible solutions. Here we want to give a general overview. The first valuable lesson we learned is to not waste time computing vanishing contributions. Many SU(2) invariants and booster functions are zero because of symmetries (triangular inequalities, odd symmetries under permutations), we avoid computing objects that are part of amplitudes that we can identify as vanishing a priori.

Symmetries also play a secondary but equally important role. The objects in the summation repeat themselves, instead of wasting time doing the same computation over and over again, we store them in the RAM and, at the end, we dump them to disk. To make it possible and not saturate the memory in a couple of iterations, we optimize the computation factoring the symmetries and consider only one representative per equivalence class. The bottleneck of the computation while dealing with SU(2) invariants is not the evaluation time but the memory needed. In an equivalence class of {6j} symbols there are approximately a hundred elements, we store only one representative per class reducing by two orders of magnitude the amount of necessary memory.

When dealing with summations with a huge amount of terms we need to be extremely careful. In our calculation, apart from the summation over auxiliary spins and intertwiners itself, these sums appear also in the evaluation of the booster functions. In particular, to deal effectively with interference effects we need to resort to arbitrary precision libraries [START_REF] Granlund | GNU multiple precision arithmetic library 6[END_REF][START_REF] Fousse | MPFR: A Multiple-Precision Binary Floating-Point Library With Correct Rounding[END_REF][START_REF] Enge | MPC: Multiple Precision Complex Library[END_REF]. How to sum a large number of finite precision floating point numbers is a well-studied problem in numerical analysis. The bound of the worst-case error on the sum of N numbers with floating point precision grows proportional to N . By using a compensated summation algorithm [START_REF] Kahan | Pracniques: Further Remarks on Reducing Truncation Errors[END_REF] we can make this error independent of N , and having it depend only on the floating-point precision of the addends. Moreover, we employ compensated summation algorithms to significantly reduce the numerical error in the addition of finite precision floating point numbers.

The full computation is divided into a large number of small tasks and we can efficiently parallelize it with OpenMP. We parallelize the computation of the booster functions and the summation over the auxiliary spins but not the evaluation of the SU(2) invariants. The computational time needed for the evaluation of the SU(2) invariants is negligible respect to all the other steps.

The sums over the six auxiliary spins are unbounded and convergent, therefore to extract a number from them we need to truncate. We introduce a homogeneous cutoff on all the spins j f ≤ l f ≤ j f + ∆ and we sum over all compatible intertwiners k e . The details of the convergence of this summation strongly depend on the value of the face spins and edge intertwiners, hindering our ability to provide a consistent estimate of the systematical error due to the truncation. We are then forced to perform an analysis a posteriori: we report some examples in Figure 8 where it is possible to estimate it between 1 and 5 percent for ∆ ≈ 5.

We mention here that we studied the possibility to develop an algorithm to guess an appropriate value for the cutoff as a function of the desired precision. However, we decided not to implement it in the current version of the library because we are lawing a way to deal with non monotonic convergences.

The error on the booster functions, due to the discretization of the integral, and the one on the SU(2) invariants are completely irrelevant with respect to the one due to the truncation. We will discuss them in more detail in the following sections.

For a single 4-simplex with fixed intertwiners, the bottleneck of the computation is the evaluation of the booster functions. If we need to compute the amplitude for all the possible intertwiners, for example, because we are contracting more than a single vertex, most of the computational time is used in performing the summation. Having in mind the goal of performing the computation of spin foam diagrams with many vertices, since the number of elements grows exponentially with it, we are studying selection rules based on the asymptotic of the booster functions or with the machine learning techniques.

SU(2) invariants

Evaluating a significant amount of SU(2) invariants can be a resources draining task. An SU(2) invariant is defined as a summation over the magnetic indices of intertwiners (and them as a sum of (3jm) symbols). Their evaluation, directly from the definition, is not convenient, as it will exhaust too many computational resources, both in time and memory. To compute an amplitude we need to evaluate a large number of invariant, requiring a technique that is both fast and memory lightweight. We can solve the time issue by reducing to one the number of summation. We can reduce every SU(2) invariant in a sum of {6j} and {9j} symbols that can be directly computed using prebuilt libraries [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 tackled this problem for the first time in [START_REF] Dona | SU(2) graph invariants, Regge actions and polytopes[END_REF] where the asymptotic of a general SU(2) invariant has been studied analytically and also numerically in the simplicial case. We employ a similar technique in sl2cfoam . Following the notation of the monography [START_REF] Yutsis | Mathematical apparatus of the theory of angular momentum[END_REF] we identify five different kinds of {15j} symbols. For example, we can write the {15j} symbols of the first kind as a sum of the product of five {6j}s symbols:

{15j} I := k 1 k 2 k 5 i 4 k 3 l 2 l 4 l 1 j 6 j 9 l 3 j 5 l 7 j 10 l 8 = x (2x + 1)(-1) i li+ji+ki+i4 k 1 l 7 x k 5 j 5 l 4 j 5 k 5 x l 8 i 4 l 6 (5)
× i 4 l 8 x k 3 j 10 j 9

j 10 k 3 x l 3 k 2 l 1 k 2 l 3 x k 1 l 7 l 2 .
To compute a single amplitude at a fixed cutoff ∆, we need to evaluate an enormous amount of {15j} symbols, and consequently of {6j} (or {9j}) symbols. We recall that all the spins l f vary in the range between the face spin j f and the cutoff j f + ∆, while the intertwiners k e assume all the values allowed by the triangular inequalities. The development of algorithms that allow for a fast and accurate computation of {6j} symbol is a fascinating and challenging topic studied by mathematicians and computer scientist, and a significant amount of libraries for the job are available. In this work, we adapt the C library 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] that can be used to compute (3jm), {6j} and {9j} symbols. We refer to their original paper for a detailed discussion of the library's performances, accuracy and memory management.

We also solved the memory usage issue by efficiently storing {6j} symbols taking into account the various symmetries of the symbols. The right data structure for this job is the hash table, a structure that efficiently maps keys to values. We associate to each {6j} symbol a key that is automatically shared by the whole class of symbols with the same value (because of symmetries). Every {6j} symbols share its value with other 143 [START_REF] Rasch | Efficient storage scheme for precalculated wigner 3j, 6j and gaunt coefficients[END_REF] symbols due to permutation symmetries. Notice that an extension of the library wigxjpf named fastwigxj provides a hash table implementation for (3jm), {6j} and {9j} symbols but, since it is designed for more general applications, it stores all possible combinations of the desired invariants with all spins smaller than a cutoff. This amount of data is way more of what we need in an ordinary EPRL four simplex computation resulting in a huge memory consumption, and consequently being an obstacle when exploring larger values of the spins. To circumvent this limitation, we implement our own hash tables based on khash provided in the library Klib [47]. We generate the keys of the hash table with the hash functions discussed in [START_REF] Rasch | Efficient storage scheme for precalculated wigner 3j, 6j and gaunt coefficients[END_REF], and we store the {6j} symbols assigning a unique six integers key to each symmetry class, saving one representative value.

We can introduce an additional simplification to improve the performance, as already used in [START_REF] Dona | SU(2) graph invariants, Regge actions and polytopes[END_REF]. We can use the freedom of choosing the intertwiners recoupling basis of the spin foam edges to rewrite the {15j} symbols as:

{15j} R := k 1 k 2 k 5 i 4 k 3 l 2 l 4 l 1 j 6 j 9 l 3 j 5 l 7 j 10 l 8 = (-1) 2(l4+l3)+2k1-(k2+k3+i4+k5) (6)
× k 1 k 3 k 2 l 1 l 2 l 3 k 1 i 4 k 5 j 6 l 4 j 5    k 2 k 3 k 1 j 10 j 9 i 4 l 7 l 8 k 4    ,
and then rewrite the {9j} symbol as a sum of three {6j} symbols:

{15j} R = (-1) 2(l4+l3)+2k1-(k2+k3+i4+k5) k 1 k 3 k 2 l 1 l 2 l 3 k 1 i 4 k 5 j 6 l 4 j 5 (7)
× x (2x + 1)(-1) 2x k 2 j 10 l 7 l 8 k 5 x k 3 j 9 l 8 j 10 x i 4 k 1 i 4 k 5 x k 2 k 3 .
Comparing with (5) we are still using five {6j} symbols but we are summing only over three of them, resulting in a relevant save of time and memory. This simplification is quite in general, (is free if the edges are not boundary ones), but we need to be careful in computations of spin foams transition amplitudes with many vertices, since the base choice on one edge affects both vertices it connects, and the split (7) needs to be done consistently in every vertex. For an explicit example of a five vertex transition amplitude where this splitting is performed see Section 4.4 of [START_REF] Dona | Infrared divergences in the EPRL-FK Spin Foam model[END_REF].

The function where we implemented the evaluation and the storing of the SU (2) invariants is called J15Symbol Hash. We sketch its working principles in the flowchart in Figure 1. The function prepares all the {6j} symbols in [START_REF] Barrett | Asymptotic analysis of Lorentzian spin foam models[END_REF] needed for the computation of a vertex amplitude given the 10 spins j f f = 1, . . . , 10 and a cutoff ∆. Following the conventions of (7) the gauge fixed edge is identified by the four spins j 5 , j 6 , j 9 and j 10 in the recoupling basis (j 5 , j 6 ) with intertwiner i 4 . After the evaluation, the hash tables are stored in a file indexed with the name of the gauge fixed edge and the cutoff ∆, with the naming convention 2j5.2j6.2j9.2j10 2i4 2∆. The factor of two multiplying the spins deals with the possible half-integer values. We decided not to use all boundary variables in the indexing to avoid files proliferation when computing amplitudes with unbounded internal faces. With our convention, all the needed symbols will be stored in the same file. In this way, we can reach a compromise between hash tables' memory consumption and time necessary to read and write the hash tables on disk. We strongly suggest selecting the gauge fixed edge as part of a boundary face to minimize the number of files created. We report a typical time and memory consumption of J15Symbol Hash in Figure 2 and As a possible upgrade for the future, we are planning to avoid the splitting of the {9j} symbol as a sum of three {6j} symbols in [START_REF] Barrett | Asymptotic analysis of Lorentzian spin foam models[END_REF] and evaluate it directly with wigxjpf with consequent hashing and storing. This upgrade hasn't been implemented yet since it requires an efficient way to assign a unique key to each symmetry class of the {9j} symbol, that we are still studying. Anyway, while this upgrade will significantly improve the speed of low spin computations, it will not be applicable in regimes with larger spins because it would use too much memory.

Booster functions

The expression for the booster functions (3) requires the calculation of the integral over the rapidity of four boost matrix elements (4) and a sum over magnetic indices. The latter is not a problem in this case since the number of terms to be summed over is relatively small (for all equal spins j is less than (2j)3 terms). However, the numerical integration is still very challenging. The integrand is a superposition of highly oscillating hypergeometric functions that, when multiplied together, interfere resulting in a localized, smooth and regular function. Moreover, for specific values, the hypergeometric functions are evaluated near one of their poles. These apparent divergences will cancel with other terms in the superposition but would require extremely high precision to control the accuracy of the result. We solved these issues by attacking them from a different angle. A direct computation is very slow, and with this motivation an alternative formula was proposed in [START_REF] Speziale | Boosting wigner's nj-symbols[END_REF] and numerically explored by [START_REF] Gozzini | Numerical study of correlations in a Lorentzian spinfoam geometry[END_REF]. Nevertheless, the situation improves significantly using a different representation for the boost matrix elements (4). This problem has been explored in [START_REF] Collet | A (simple) expression of the unitary-irreducible SL(2,C) representations as a finite sum of exponentials[END_REF] where each boost matrix elements (4) is rewritten as a finite sum of exponentials with complex coefficients. For simple irreducible representations and γ = 0, the result is given by: where Y are complex coefficients 3 . In the same work they provide a C++ implementation of the booster functions based on this formula. We integrated this code in sl2cfoam porting it in C, we added the evaluation of the (4jm) symbols in (3) with wigxjpf and we provided an hash table data structure similar to the one we introduced for the SU(2) invariants in the previous section. The evaluation of these functions requires arbitrary precision floating point numbers, implemented with MPFR and MPC. These tools allow us to compute the booster functions in an accurate way and to test some of theirs properties.

d (γj,j) jlp (r) = 1 (e r -e -r ) j+l+1 (8)
The integral is performed using the trapezoidal rule with 3000 points as a default number. We explored this choice in Figure 4. The convergence looks very fast in the number of points; we observe a slight worsening with the increase of the l f . The number of integration points can be increased arbitrarily to reduce even further the admissible error on the single booster function. The sampling of the integration interval is parallelized with OpenMP. The library sl2cfoam implement the evaluation of all the booster functions needed in the computation of a vertex amplitude with the function B4 Hash. We illustrate its structure in the diagram in Figure 5.

The values of the booster functions are stored in hash tables and then in a file indexed with the name of the gauge fixed edge and the cutoff ∆, using the same naming conventions as the {6j} symbols 2j5.2j6.2j9.2j10 2∆. Furthermore, we divide the files into subfolders depending on the Immirzi parameter γ to avoid mixing of incompatible evaluations in case one is interested in working with different values for γ. We store in each file the booster functions with all the possible combinations of the six l f spins appearing in a vertex amplitude and every compatible intertwiner. We report a typical time consumption of B4 Hash in Figure 6 and the scaling of the time as a function of the number of cores used in the computation in Figure 7. The memory's use of booster functions is negligible compared to {6j} symbols' one.

For completeness, we provide also a function, called B4Function, that computes a single booster, in case one is interested in them outside the context of vertex amplitudes.

When computing the booster functions, one notices that, for a fixed set of spins j f and l f , not all combinations of intertwiners contributes with values of similar order of magnitude. An example of this can be seen in [START_REF] Speziale | Boosting wigner's nj-symbols[END_REF] and in [START_REF] Puchta | Asymptotic of Lorentzian Polyhedra Propagator[END_REF], for the case l f = j f (called simplified model), where the main contribution is given by the same intertwiners i e = k e . The generalization of this result for any set of l f is unfortunately still unknown. This result, or an equivalent method to localize the dominant contribution from the booster functions, will immediately translate into a significant save of time in the evaluation of the amplitude, especially in the case of multiple vertices.

We also plan to explore the possibility of completely changing the way in which we compute the booster functions and use their expression in terms of the Clebsch-Gordan coefficients for SL(2, C) introduced in [START_REF] Speziale | Boosting wigner's nj-symbols[END_REF]. Using the Clebsch-Gordan formulation of the booster functions, we trade the integrals over the rapidity and the sum over the magnetic indices, with a single integral (and a finite sum) over a unitary SL(2, C) virtual irrep. To do so, one has to compute Clebsch-Gordan coefficients fast and with high accuracy. The speed of the functions we included in sl2cfoam that do this task is too low to compete with the actual implementation of the booster functions. For this reason, we have not implemented it yet. The possibility of using recursion relations fro the Clebsch-Gordan coefficients also makes this option very appealing. A first attempt was made in [START_REF] Gozzini | Numerical study of correlations in a Lorentzian spinfoam geometry[END_REF], we think this is a very promising route to improve sl2cfoam .

B4 Hash

Four simplex amplitude

In the previous two sections we illustrated the evaluation of the constituents of the vertex amplitude [START_REF] Etera | A New spinfoam vertex for quantum gravity[END_REF]. In this section we will discuss how all the contributions are summed together. The summation is performed over six auxiliary spins l f and four k e intertwiners, whose ranges depends on the l f but are generally bounded. The main problem in this step is the extension of these summations. They are, generally, unbounded and numerous. As we already stressed multiple times, the full amplitude is well defined [START_REF] Engle | Regularization and finiteness of the Lorentzian LQG vertices[END_REF], implying that the summations over the l f are convergent, thus being unbounded is not a problem. Nevertheless, to be able to do a numerical evaluation of the sum we need to introduce a homogeneous cut-off ∆. In the following, we will denote as the N th shell the set of factors contributing to the sum with maximum spin difference Max {l fj f } = N . To be concrete we studied the convergence of the amplitude as a function of the cutoff, generally, it is quite rapid. We report some explicit examples in some particular cases in Figure 8. It is not possible to have a unique prescription to set the optimal ∆, since the convergence of the amplitude depends on the details of the data like the face spins j f and the Immirzi parameter. We explored other configurations and the convergence is qualitatively the same. However, the ratio between the amplitude value with a different number of shells varies. In the library sl2cfoam the function FourSimplex performs the summation of the 4-simplex vertex amplitude. The algorithm at the base of the implementation is described in Figure 12. This function is the main result of this work and the big star of the library. Given the ten face spins, the five edge intertwiners, the cutoff ∆ and the Immirzi parameter, FourSimplex uses the functions J15Symbol Hash and B4 Hash described in the previous sections to precompute all the needed SU(2) invariants and booster functions and then proceed to the summations. This last step needs particular care. Given the vast number of factors to be summed over we need to employ a compensated summation algorithm to reduce the numerical error caused by the sum of finite precision floating point numbers. Furthermore, to minimize the error, we convert and multiply the components of the factors to arbitrary precision floating point numbers (using the library MPFR). Parallelization is needed at this step but, while we have implemented it for the summation over the face spins l f , it was counterproductive to do it also for the four intertwiner sums k e because of triangular inequalities: when we fix a set of l f we set the bounds over the k e and the remaining four sums depend on them. In Figure 9 we show a parallelization test at fixed workload while increasing the number of cores.

Figure 9: Parallelization at fixed workload: Scaling of the time required by the computation of a 4-simplex configuration with j f = 1 and ∆ = 11 while increasing the number of cores. The power law a × j b is given by a = 109s and b = -0.76. Since the parallelization acts only on the l f summations and not on the ke ones the code performs better, in this regard, when a large ∆ is considered.

The data are saved using the same convention as the boosters functions, we index them as 2j5.2j6.2j9.2j10 2∆. We stress that selection rules, to find the important terms within the sums, would be an essential step forward. At the present stage, we can set an empirical cut-off over boosters to reduce the number of terms to be summed over but it would be exciting to study, with new methods such as machine learning's techniques, selection rules for the vertex amplitude. Evaluation time of a single factor in a computation with all boundary spins and intertwiners equal to λ. Each factor takes ∼ 10 -5 s. Memory data access using khash works at ∼ 10 -7 s: we have to retrieve data for four B4, two "fixed" {6j} and for a sum of three {6j}, converting data to MPFR variables and then assembly the amplitude.This leads us to ∼ 10 -5 s of evaluation time for a single factor of the l f , ke sum. If we would be able in the future to directly store {9j} symbols we could save one order of magnitude in time. The study of the semiclassical behavior of a spin foam transition amplitude relies on the relation between polyhedra and SU(2) invariants [START_REF] Barrett | Asymptotic analysis of the eprl four-simplex amplitude[END_REF][START_REF] Dona | SU(2) graph invariants, Regge actions and polytopes[END_REF] and coherent intertwiners [START_REF] Etera | A New spinfoam vertex for quantum gravity[END_REF]. Thus, in sl2cfoam , we provide a tool to compute Livine-Speziale coherent intertwiners as arbitrary precision floating point using MPFR. We refer to [START_REF] Dona | SU(2) graph invariants, Regge actions and polytopes[END_REF] for the conventions and background material used in the present work, while we report here the definitions to fix the notation. We recall that an n-valent coherent intertwiner [START_REF] Etera | A New spinfoam vertex for quantum gravity[END_REF] is defined by the group averaging:

FourSimplex

| |{j a , n a } := dg ⊗ n a=1 g |j a , n a ∈ Inv ⊗ n a=1 V (ja) , (9)
where V ja are SU(2) irreducible representation, Inv ⊗ n a=1 V (ja) is a singlet vector space and the vectors n a can be parametrized as n a := (sin Θ cos Φ, sin Θ sin Φ, cos Θ). The library provides a function to compute them for n = 4. In this case the decomposition of a coherent intertwiner, in the recoupling channel i with outgoing links, reads

c i ( n a ) j a , i| j a , n a = ma j 1 j 2 j 3 j 4 m 1 m 2 m 3 m 4 (i) 4 a=1 D ja ma,ja ( n a ) = ñ1 ñ2 ñ3 ñ4 j 1 j 2 j 3 j 4 i , (10)
where D j m,j ( n) = D j m,j (φ, Θ, -φ) is the SU(2) Wigner matrix parametrized as:

D j m,j (φ, Θ, -φ) = e -imφ d j m,j (Θ)e ijφ ,
with d j m,j as a small Wigner matrix [START_REF] Varshalovich | Quantum Theory of Angular Momentum[END_REF]. CoherentStateI computes a coherent state coefficient c i ( n a ) for a given intertwiner i. It takes as input the four spins, one intertwiner, the four set of angles (Θ a , φ a ) (parametrizing the four normals n a ) and also four signs related to the orientation of the links. For instance, to compute the coherent intertwiner in [START_REF] Bianchi | LQG propagator from the new spin foams[END_REF] one needs to set:

CoherentStateI (2i, 2j 1 , 2j 2 , 2j 3 , 2j 4 , φ 1 , Θ 1 , φ 2 , Θ 2 , φ 3 , Θ 3 , φ 4 , Θ 4 , +1, +1, +1, +1) .
The last four arguments represent the "outgoing" orientation of the four links in the same order the spins are inputted. To represent "ingoing" links is sufficient to flip the sign to -1 of the corresponding entries.

SL(2, C) symbols

In this section we review the basic formulas of Lorentzian recoupling theory and introduce the code of sl2cfoam to compute them. We will follow closely the notation in [START_REF] Speziale | Boosting wigner's nj-symbols[END_REF]. The SL(2, C) Clebsch-Gordan coefficients can be factorized in an SU(2) one and a boost one:

C ρ3k3j3m3 ρ1k1j1m1ρ2k2j2m2 = χ(ρ 1 , ρ 2 , ρ 3 , k 1 , k 2 , k 3 ; j 1 , j 2 , j 3 ) C j3m3 j1m1j2m2 . (11)
We will use a simplified notation indicating with χ(j 1 , j 2 , j 3 ) = χ(ρ 1 , ρ 2 , ρ 3 , k 1 , k 2 , k 3 ; j 1 , j 2 , j 3 ). Denoting with

J = i j i , K = i k i , P = i ρ i , (12)
we can write the boost part of the SL(2, C) Clebsch-Gordan coefficients as a finite sum of ratios of complex gamma functions as:

χ(j 1 , j 2 , j 3 ) = (-1) K+J 2 4 √ 2π x (ρ i , k i )Γ( 1-iP +K 2 )Γ( 1-iP -K 2 ) d j1 d j2 d j3 κ(ρ i , k i ; j i ), (13)
with

κ(ρ i , k i ; j i ) = (-1) -k1-k2 e -i Φ ρ 1 j 1 +Φ ρ 2 j 2 -Φ ρ 3 j 3 (-1) j1-j2+j3 d j3 (j 1 -k 1 )!(j 2 + k 2 )! (j 1 + k 1 )!(j 2 -k 2 )! 1/2 (14) × min{j1,k3+j2} n=-j1 (j 1 -n)!(j 2 + k 3 -n)! (j 1 + n)!(j 2 -k 3 + n)! 1/2 j 1 j 2 j 3 n k 3 -n -k 3 × j1 s1=max{k1,n} j2 s2=max{-k2,n-k3} (-1) s1+s2-k1+k2 (j 1 + s 1 )!(j 2 + s 2 )! (j 1 -s 1 )!(s 1 -k 1 )!(s 1 -n)!(j 2 -s 2 )!(s 2 + k 2 )!(k 3 -n + s 2 )! × Γ( 1-i(ρ1-ρ2-ρ3)-K+2s1 2 )Γ( 1+i(ρ1-ρ2+ρ3)+K+2s2 2 )Γ( 1-i(ρ1+ρ2-ρ3)-k1+k2+k3-2n+2s1+2s2 2 ) Γ(1 -iρ 1 + s 1 )Γ(1 -iρ 2 + s 2 )Γ(1 + iρ 3 + s 1 + s 2 )Γ( 1-iP -k1+k2+k3-2n 2
) ,

where the additional phase

x (ρ i , k i ) = Γ 1+iP -K) 2 |Γ 1+iP -K) 2 | (15) × Γ 1-i(-ρ1+ρ2+ρ3)-k1+k2+k3 2 |Γ 1-i(-ρ1+ρ2+ρ3)-k1+k2+k3 2 | Γ 1-i(ρ1-ρ2+ρ3)+k1-k2+k3 2 |Γ 1-i(ρ1-ρ2+ρ3)+k1-k2+k3 2 | Γ 1-i(-ρ1-ρ2+ρ3)-k1-k2+k3 2 |Γ 1-i(-ρ1-ρ2+ρ3)-k1-k2+k3 2 | ,
was introduced in [START_REF] Speziale | Boosting wigner's nj-symbols[END_REF] to make the Clebsch-Gordan coefficients, and consequently the booster functions, always real. This phase is different from the convention used in [START_REF] Ruhl | The Lorentz group and harmonic analysis[END_REF].

We can write a SL(2, C) version of the (3jm) symbol. It reads

ρ 1 , k 1 ρ 2 , k 2 ρ 3 , k 3 j 1 j 2 j 3 = (-1) j1-j2+j3 d j3 χ(j 1 , j 2 , j 3 ). (16)
We provide a function, wigner 3j sl2c, to compute it. Similarly we can also define an analog of the SU(2) {6j} symbol for SL(2, C) as:

{6ρ, 6k} = ji≥ki {6j i } ρ 1 , k 1 ρ 2 , k 2 ρ 3 , k 3 j 1 j 2 j 3 ρ 1 , k 1 ρ 5 , k 5 ρ 6 , k 6 j 1 j 5 j 6 (17)
× ρ 4 , k 4 ρ 2 , k 2 ρ 6 , k 6 j 4 j 2 j 6 ρ 4 , k 4 ρ 5 , k 5 ρ 3 , k 3 j 4 j 5 j 3 ,
To perform a numerical evaluation of this symbol we artifically cutoff the summation to j i ≤ ∆. We provide wigner 6j sl2c to compute it. Since the κ in ( 14) is given by a sum of ratio of complex gamma functions, an arbitrary precision computation is needed. We implemented it using MPFR variables. We tested the commands of this section up to spins ≈ 50, we leave the improvement of this upper limit to future works.

What comes next?

In this paper, we described the library sl2cfoam that we developed to compute the spin foam vertex amplitude defined by the EPRL model. Many technical improvements are necessary to unlock the full potential of this code. In analogy to what we do for the {6j} symbols, we need to implement hash tables also for {9j} symbols. This would speed up the computation of almost an order of magnitude, trading at low spins memory usage with faster access to data. Moreover, we want to explore the possibility of computing booster functions with the Kerimov formula, as explained in [START_REF] Speziale | Boosting wigner's nj-symbols[END_REF], for which recursion relations are available, simplifying the summation over the auxiliary spins l f . To perform the computation of transition amplitudes with many vertices, it is crucial to improve our parallelization process, using MPI and POSIX threads enabling many servers calculations. A similar improvement will be reached by finding selection rules to restrict the summation to only the important terms. Our lack of analytical insight into the booster functions will force us to use new techniques like for example machine learning or Monte Carlo techniques. Finally, we want to search for a dynamical way of dealing with the cutoff ∆ to be able to perform computations with an error set a priori. We plan to explore these fascinating possibilities soon. For a quick reference, we report in Table 1 the timing of the various functions in sl2cfoam . We firmly believe that this code could profoundly change the way we understand the dynamics of loop quantum gravity. One of the most successful computations in the theory proved the connection in the semiclassical limit with discrete General Relativity [START_REF] Barrett | Asymptotic analysis of the eprl four-simplex amplitude[END_REF][START_REF] Barrett | Asymptotic analysis of Lorentzian spin foam models[END_REF]. They showed that a single EPRL vertex amplitude contains the exponential of Regge action for a Lorentzian 4-simplex in the large spin limit. While this calculation is of dramatic importance for the theory, its reign of validity is limited to a single simplex. We still lack a definite result for an extended triangulation made of many simplices glued together. This limitation made the community questions if the simplicity constraints, responsible of constraining the original topological BF theory to General Relativity, are enforced correctly [START_REF] Conrady | On the semiclassical limit of 4d spin foam models[END_REF][START_REF] Bonzom | Spin foam models for quantum gravity from lattice path integrals[END_REF][START_REF] Hellmann | Holonomy spin foam models: Asymptotic geometry of the partition function[END_REF]. This translate to an inability of the model in describing amplitudes compatible with curved spacetime. Preliminary results, limited to the Euclidean version of the model, show that curved geometries are indeed allowed in the asymptotic limit of the amplitude. Using sl2cfoam we will investigate the large spin limit of the ∆ 3 complex, dual to three 4-simplices sharing a bulk triangle, the most straightforward triangulation with an internal face.

j f = i 4 = 7 j f = i 4 = 15 j f = 1, i 4 = 1, ∆ = 5 j f = 1, i 4 = 1, ∆ =
We are performing a consistency check of our code by computing the homogeneous scaling of the Lorentzian single vertex amplitude with Euclidean boundary data. To do so, we calculated the vertex amplitude with spins, both integers and half-integers up to order 10, contracted with coherent intertwiner states representing the boundary of a regular Euclidean four simplex. We sum two shells per amplitude with an estimated error of approximately the 20%. We find an amazing agreement with the analytic formula of [START_REF] Barrett | Asymptotic analysis of the eprl four-simplex amplitude[END_REF].

We also plan to explore numerically the divergence of the theory recently estimated in [START_REF] Riello | Self-energy of the Lorentzian Engle-Pereira-Rovelli-Livine and Freidel-Krasnov model of quantum gravity[END_REF][START_REF] Dona | Infrared divergences in the EPRL-FK Spin Foam model[END_REF] by computing the "self-energy" diagram. A large set of power-law divergences is still compatible with both estimates, and the tensorial structure of the diagram is still unknown. To confirm the first and compute the second will allow us to understand how to correctly regularize the theory without introducing a cosmological constant following a prescription similar to the one studied in [START_REF] Freidel | Diffeomorphisms and spin foam models[END_REF].

The renormalization group flow of the EPRL vertex amplitude is of crucial importance for the future of the theory. Finding a fixed point in the flow and identifying phase transitions would allow us to understand if the diffeomorphism invariance, broken by the discretization, is recovered and what mechanism will restore it. An interesting insight on this topic is provided in [START_REF] Dittrich | The continuum limit of loop quantum gravity -a framework for solving the theory[END_REF][START_REF] Dittrich | Coarse graining flow of spin foam intertwiners[END_REF][START_REF] Delcamp | Towards a phase diagram for spin foams[END_REF] and in [START_REF] Bahr | On background-independent renormalization of spin foam models[END_REF][START_REF] Bahr | Numerical evidence for a phase transition in 4d spin foam quantum gravity[END_REF][START_REF] Bahr | Hypercuboidal renormalization in spin foam quantum gravity[END_REF][START_REF] Bahr | Towards a Cosmological subsector of Spin Foam Quantum Gravity[END_REF][START_REF] Bahr | Renormalization in symmetry restricted spin foam models with curvature[END_REF]. In this last series of works they perform a numerical investigation of the flow of the Euclidean EPRL model and on a hypercubic lattice, with a large number of vertices. The calculations are possible thanks to some strong but vital approximations: they work with the Euclidean version of the theory, they restrict to a specific form of the intertwiners, and instead of using the full vertex amplitude they replace it with its asymptotic formula.

At the present sl2cfoam cannot directly attack analog computation, the number of vertices needed for a lattice computation is still out of our reach, but we do not exclude to be able to do it in a near future after all the proposed upgrade are implemented and we have access to more processing power (for the moment we are limited to a 32 cores machine).

Nevertheless, if we restrict to a purely quantum region characterized with relatively small spins, probably up to order 10, we can now compute any transition amplitude with a few vertices. We would like to invite the community to look for interesting computations that could be performed within those limitations, and to develop a common framework in which the theory can be used to answer physically meaningful questions.

A SU(2) Symbols

We use the definition of the Wigner's (3jm) symbol reported in [START_REF] Varshalovich | Quantum Theory of Angular Momentum[END_REF] with the following orthogonality properties m1,m2

j 1 j 2 j 3 m 1 m 2 m 3 j 1 j 2 j 3 m 1 m 2 n 3 = δ j3l3 δ m3n3 2j 3 + 1 , (18)
implying they are normalized to one and that, if triangular inequalities are not satisfied, they vanish. We can define the (4jm) symbol in the following way, as the contraction of two (3jm) symbol using an intertwiner i

with the orthogonality relations m1,m2,m3

j 1 j 2 j 3 j 4 m 1 m 2 m 3 m 4 (i1) j 1 j 2 j 3 j 4 m 1 m 2 m 3 m 4 (i2) = δ i1i2 2i 1 + 1 δ j4l4 δ m4n4 2j 4 + 1 , (20)
where

δi 1 i 2
2i1+1 is the normalization factor.

Figure 1 :

 1 Figure 1: J15Symbol Hash is the main function to store in an efficient way {6j} symbols.

Figure 2 :

 2 Figure 2: Evaluation time of the function J15Symbol Hash. Left panel: Scaling of the configuration with all thespins and the gauge fixed intertwiner equal j f = i4 = λ and a vanishing cutoff ∆ = 0, we compute the {15j} symbols for all the values of the intertwiners ke. We see a power law trend a × j b with a = 2.9 × 10 -6 s and b = 5.3. Right panel: Scaling in the cutoff of the configuration with all the spins and the gauge fixed intertwiner equal to j f = i4 = 1. We see a power law trend a × j b with a = 4 × 10 -4 s and b = 5.0.

Figure 3 :

 3 Figure 3: Memory usage of the function J15Symbol Hash. Left panel: Scaling of the configuration with all the spins and the gauge fixed intertwiner equal j f = i4 = λ and a vanishing cutoff ∆ = 0, we compute the {15j} symbols for all the values of the intertwiners ke. We see a power law trend a × j b with a = 1.8 × 10 -4 Mb and b = 5.1. Right panel: Scaling in the cutoff of the configuration with all the spins and the gauge fixed intertwiner equal to j f = i4 = 1. We see a power law trend a × j b with a = 1.8 × 10 -3 Mb and b = 3.7.

Y

  γj,j;j,l;p m e (j+l-|j-p|-2m-iγj)r + (-1)

Figure 4 :

 4 Figure 4: Dependence of the booster function on the number of points in the rapidity integral. We plot the relative difference ǫ with respect to the value with 10 4 points. All three booster functions reported here have all equal spins j f = 10 and intertwiner i = 0. The four spins l f and intertwiner k are respectively: l f = 11, 12, 13, 10 and k = 1, l f = 20, 22, 23, 20 and k = 2, l f = 30, 32, 33, 30 and k = 2.

Figure 5 :

 5 Figure 5: B4 Hash compute and store efficiently the booster functions.

Figure 6 :

 6 Figure 6: Evaluation time of the function B4 Hash. Left panel: Scaling of the configuration with all the spins andthe gauge fixed intertwiner equal j f = i4 = λ and a vanishing cutoff ∆ = 0. We compute all the booster functions for all the values of the intertwiners ke. We see a power law trend a × j b with a = 0.7s and b = 2.5. Right panel: Scaling in the cutoff of the configuration with all the spins and the gauge fixed intertwiner equal to j f = i4 = 1. We see a power law trend a × j b with a = 2.7s and b = 4.3.

Figure 7 :

 7 Figure 7: Parallelization at fixed workload: scaling of the time required by the computation of booster functions for a 4-simplex with j f = 10 and ∆ = 0 while increasing the number of cores. The power law a × j b reads a = 209s and b = -0.7.

Figure 8 :

 8 Figure 8: Convergence of the vertex amplitude. Left panel: Convergence of the vertex amplitude while increasing the cutoff ∆ for a configuration with all j f = 1 and three different set of boundary intertwiners. Right panel: The relative error ǫ between the shell ∆ and the shell ∆ = 12, that has been taken as the correct value of the amplitude. The error for the simplified model (∆ = 0) is ∼ 50% while from ∆ = 5 it decreases to 1%-5%. The convergence of the vertex amplitude depends on many factors including the chosen boundary configuration and the Immirzi parameter.

Figure 10 :

 10 Figure 10: Evaluation time of the function FourSimplex. Left panel: Scaling of the configuration with all the spins and the gauge fixed intertwiner equal j f = i4 = λ and ∆ = 0. We see a power law trend a × j b with a = 2.0 × 10 -6 s and b = 4.8. Right panel: Scaling in the cutoff of the configuration with all the spins and the gauge fixed intertwiner equal to j f = i4 = 1. We see a power law trend a × j b with a = 4 × 10 -4 s and b = 5.2.

Figure 11 :

 11 Figure 11: Left panel: We show the number of factors # f act to be cycled over to a single 4-simplex with boundary spins, intertwiners equal to λ, and ∆ = 0. The power law is given by a = 0.2 and b = 4.3 Right panel:Evaluation time of a single factor in a computation with all boundary spins and intertwiners equal to λ. Each factor takes ∼ 10 -5 s. Memory data access using khash works at ∼ 10 -7 s: we have to retrieve data for four B4, two "fixed" {6j} and for a sum of three {6j}, converting data to MPFR variables and then assembly the amplitude.This leads us to ∼ 10 -5 s of evaluation time for a single factor of the l f , ke sum. If we would be able in the future to directly store {9j} symbols we could save one order of magnitude in time.

  3.

						input: ten spins,
				J15Symbol Hash homogeneous cut off ∆
	yes	Have you the table with the required ∆ ?	yes	Have you f6jg symbol hash tables?	no	Initialize a new hash table
							cycle all possible
							combinations of l f and ke
		no
			Load the table with	cycle all possible combinations of l f and ke	Fix a set of six li spins
		the highest ∆ available			and four ke intertwiners
							f6jg key conversion
	End		Does the key already exist?		no		compute and store the f6jg value Create a new key,
	Store the hash table			Is this the last set?

  FourSimplex is the main function to compute and store EPRL vertex amplitudes

	Return value End Store and return value no Have you the table with the required ∆ and the required value? Figure 12: 4 Utilities yes 4.1 Coherent states	yes yes	Fix a set of six l f spins and four ke intertwiners Have you 4-simplex hash tables? Is this the last set? no	homogeneous cut off ∆ input: ten spins, Initialize a new hash table no J15Symbol Hash Immirzi parameter five intertwiners B4 Hash parallelization on l f Compute value and sum Load f6jg and boosters tables start summation conversion in MPFR variables compensated summation

Table 1 :

 1 Summary of the computational time of the functions J15Symbol Hash, B4 Hash, and FourSimplex. We report some examples of computational time for configuration with all spins equal to 7 and 15 in the first two columns and the parameters of the power law behavior aj b . In the last two columns we report some examples of computational time for configuration with all spins equal to 1 as a function of the cutoff and the parameters of the power law behavior a∆ b . All the computations are performed with a single processor.

	10

The variables of Regge calculus are the simplex lengths (encoding the intrinsic curvature). The extrinsic curvature is encoded in the dihedral angles hinging on the triangles of the triangulation.

that are effectively magnetic indices respect the group SL(2, C) and are originated by the splitting of two representation matrices

We refer to[START_REF] Collet | A (simple) expression of the unitary-irreducible SL(2,C) representations as a finite sum of exponentials[END_REF] for their full expression.

Acknowledgments

This work was supported in part by the NSF grants PHY-1505411, PHY-1806356 and the Eberly research funds of Penn State. We want to thank François Collet and Francesco Gozzini for providing a basic code for the booster function and the complex gamma function with arbitrary precision floating points numbers and especially for many interesting and helpful discussions. Thanks to NORDITA and to the organizers of the workshop "Quantum Gravity on the Computer", a very inspiring and useful week in Stockholm. Many thanks to Simone Speziale for fruitful discussions and encouragement and we are very happy to welcome Emilio in this puzzling world.