How-To compute EPRL spin foam amplitudes

Spin foam theory is a concrete framework for quantum gravity where numerical calculations of transition amplitudes are possible. Recently, the field became very active, but the entry barrier is steep, mainly because of its unusual language and notions scattered around the literature. This paper is a pedagogical guide to spin foam transition amplitude calculations. We show how to write an EPRL-FK transition amplitude, from the definition of the 2-complex to its numerical implementation using \texttt{sl2cfoam-next}. We guide the reader using an explicit example balancing mathematical rigor with a practical approach. We discuss the advantages and disadvantages of this approach and provide a novel look at a recently proposed approximation scheme.


Introduction
Spin foam theory provides a background-independent, Lorentz covariant path integral for general relativity. Spin foams provide dynamics to Loop Quantum Gravity, defining transition amplitudes between spin network states. A triangulation discretizes the space-time manifold, and its 2-complex regularizes the partition function.
The EPRL-FK model [1,2] (we will refer to it as just EPRL for brevity) is the state-of-the-art spin foam model. There is a large consensus in the community [3,4,5,6,7] that the classical continuum theory can be recovered with a double limit of finer discretization and vanishing . This observation is supported by the emergence of Regge geometries and the Regge action in the asymptotics of the 4-simplex vertex amplitude for large quantum numbers [8,9] and the recent study of many vertices transition amplitudes.
The amount of calculations possible within the models recently grew considerably. It was possible because of a paradigm shift in the field. It evolved from a theoretical framework to circumvent the difficulties in imposing the Hamiltonian constraint in the canonical approach [10] into a concrete tool where numerical calculations of transition amplitudes are possible [11,12,13,14,15,16,17].
With increased interest in the field, its entry barrier also increased vastly. Getting into spin foam is very difficult for a student or a researcher from a different field. There are plenty of reviews [18,19] and excellent books [20] to study and learn the basic theory. On the other end of the spectrum, we have plenty of advanced frontline papers that explore the connection of spin foam with GR [3,21,6] or possible phenomenological implications [22,23,24,25].
We noticed a hole in the literature. There are no papers that give you all the tools needed to complete a spin foam calculation, from its conception to the number. With this paper, we guide the reader through the calculation of an EPRL transition amplitude in a pedagogical manner. We use an explicit example to help them not feel disoriented dealing with abstract concepts. We hope that this paper can fill that hole and open spin foam calculation to a new generation of students and researchers.
To read this paper, advanced background knowledge on spin foam is not necessary. However, a basic understanding of the topic is helpful. We think of this work as a guide for making spin foam calculations. We refer to targeted reviews of the EPRL model [18,20,19] for a comprehensive discussion of its definition, motivation, and physical significance.
We start with a brief review of the construction of the spin foam theory and the definition of the EPRL model in Section 2. In the rest of the paper, we show the reader how to compute a spin foam transition amplitude associated with a triangulation of the space-time manifold. We identify five necessary steps, each illustrated in a different section.
In Section 3, we describe how to build the 2-complex from the triangulation. It is crucial in regularizing the gravitational path integral and writing a finite transition amplitude.
Step 2. Write the EPRL spin foam amplitude.
In Section 4, we give the prescription to write the transition amplitude associated with a 2complex, and we introduce a very convenient graphical method to represent the amplitude. For the calculation, we resort to a divide-and-conquer strategy.
Step 3. Divide the EPRL transition amplitude into vertex contributions.
In Section 5, we show how to divide any transition amplitude into vertex amplitudes.
Step 4. Compute the EPRL vertex amplitudes. In Section 6, we discuss the calculation of the vertex amplitude in terms of SU (2) invariants and booster functions.
Step 5. Use sl2cfoam-next to compute a number. We perform the numerical evaluation of the amplitude in Section 7 using the numerical library sl2cfoam-next and discuss the necessary approximations. In this section, we also discuss and improve the extrapolation scheme discussed in [13] as a tentative to lift, at least part of, the approximation used to calculate the amplitude.
We complement our discussion with an explicit example. We compute the EPRL transition amplitude based on the triangulation ∆ 4 . It was considered first in [26] in Lorentzian the spin foams. It is 2-complex that at the same time not trivial (with more than one vertex), simple (with four vertices and some symmetry) but rich enough (with one bulk face) to require a certain degree of optimization to compute the associated amplitude. Moreover, in [26] coherent boundary data corresponding to a Lorentzian geometry was provided, allowing semiclassical calculation with some ease that we leave to future work.

Overview of the EPRL model
Spin foam theory is a promising approach to quantize gravity. The goal is to define a path integral for general relativity in a non-perturbative and background-independent way. The spin foam partition function assigns transition amplitudes between spin network states, a basis of the Loop Quantum Gravity kinematical Hilbert space. For this reason, spin foam theory gives a dynamic to Loop Quantum Gravity, and it is often referred to as Covariant Loop Quantum Gravity [20].
In the Plebanski formulation of general relativity [27], we formulate gravity as a topological BF theory with constraints [28]. The variables of the BF theory are a 2-form B conjugated to a connection ω (with curvature F ) 1 .
General relativity is not topological, the simplictiy constraints reduce the B-field in BF theory to a γ-simple 2-form B = ⋆e ∧ e + 1 γ e ∧ e , reducing the action to the familiar Holst action [29]. The path integral of spin foam theory is regularized on a triangulation, more precisely its 2-complex, to truncate the degrees of freedom. We discretize and quantize the topological theory first. The B-fields are assigned to faces of the 2-complex, triangles, and encode their geometry. The connection is regularized by considering only its holonomy g responsible for the parallel transport along the (half-)edges of the 2-complex (from one tetrahedra to another). The topological theory partition function consists of a collection of delta functions imposing flatness of each face of the 2-complex.
The partition function of the EPRL model is derived enforcing the simplicity constraints at the quantum level to reduce the topological theory to gravity. On a simplicial triangulation, we have a linear version of the simplicity constraints: we require the proportionality between the boost and rotation generators of SL(2, C) K = γ L at the boundary of any 4-simplex (vertex of the 2-complex). The generalization to arbitrary tessellation is possible [30] but requires complications beyond this work's scope. Therefore, we limit ourselves to 4-simplices.
The key ingredient of the EPRL model is the Y γ map. It embeds the spin j SU (2) representation into the lowest spin sector of the unitary irreducible representations in the principal series of SL(2, C) labeled by ρ, k = γj, j.
We expand the BF theory partition function in terms of matrix elements of the holonomies in irreducible rerpesentations of SL(2, C) D ρ,k jmln (g). See Appendix B and references therein for more details. The EPRL model prescription enforces the Y γ map at every vertex of the 2-complex restricts the irreducible representations to γ-simple ones D γj,j jmjn (g) [1]. If the 2-complex has a boundary, the spin foam partition function maps states from the Loop quantum Gravity kinematical Hilbert space (identified with the boundary space of the spin foam with the Y γ map) into the complex numbers (quantum transition amplitudes between these states).
The EPRL spin foam partition function is given as a state sum over SU (2) spins j f on the faces and intertwiners i e on the edges of the 2-complex: 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 [32]. We will not give an explicit form of the amplitude for an arbitrary 2-complex. They can be found in many references [1,18,20] if the reader is interested. Instead, we opt for a constructive approach. In Section 4, we guide the reader through a set of rules to write a general EPRL transition amplitude. In Section 5, we divide the transition amplitude in vertex amplitudes. In Section 6, we discuss the explicit form of the vertex amplitude and its form best suited for numerical calculations.

How-To draw the 2-complex
The spin foam partition function is regularized on the 2-complex of a triangulation of the space-time manifold. Given a triangulation, we can build its 2-complex associating a vertex to each 4-simplex.
Each 4-simplex shares a tetrahedron with an adjacent 4-simplex. We associate to each tetrahedron of the triangulation an edge of the 2-complex. An edge connects two adjacent vertices. Each triangle in a 4-simplex is shared by two tetrahedra, which are generally shared with other 4-simplices. In the whole triangulation, a triangle can be shared by any number of tetrahedra and 4-simplices.
We associate to each triangle of the triangulation a face of the 2-complex. A face can contain any number of vertices and all the edges connecting them. We also assign an orientation to the faces of the 2-complex. This choice is needed for a well-defined notion of parallel transport (to identify the source and target of the holonomy uniquely).
Since two vertices share each edge for each of them, we can introduce two half-edges, one associated with each vertex. We can still picture them as dual to the tetrahedron but "seen" in the 4-simplex it belongs to. In each vertex in a given face there are two half-edges. This is sometimes referred to as a wedge. The orientation of the face allows us to identify one half-edge as the source tetrahedron (reference frame) and the other as the target tetrahedron (reference frame) of parallel transport along the face from the first one to the last one. If a boundary is present, the edges intersected by the boundary are severed in half, leaving only one half-edge in the skeleton.
We summarize the nomenclature introduced in this section in the following:

An example: the ∆ 4 triangulation
The triangulation is formed by four 4-simplices, all sharing a triangle. The triangulation has seven points, nineteen segments, twenty-five triangles, sixteen tetrahedra (twelve in the boundary and four in the bulk), and four 4-simplices. We label the points with numbers from 1 to 7, segments with couples of different numbers (points), triangles with triples of distinct numbers (the shared triangle is 123), tetrahedra with a quadruple of distinct numbers, and 4-simplices with five distinct numbers. See Figure 1 for a pictorial representation of the triangulation.  edges, each associated with a tetrahedron shared among two 4-simplices. There are also three external edges for each vertex. Each edge belongs to 4 faces, and each face is associated with a triangle of the ∆ 4 triangulation. All triangles but one belong to the boundary of the triangulation. Therefore all faces but one of the 2-complex are boundary faces. The bulk face is associated with the triangle shared by all the 4-simplices. Thus it is crossing all four vertices. We label the 2-complex in the same way of the triangulation, see Figure 2 for a representation.

How-To write the EPRL spin foam amplitude
For each wedge we write a γ-simple unitary irreducible representation in the principal series of SL(2, C) (see Appendix B and reviews [31] and references therein for more mathematical details).
where j ∈ N/2 is a spin, γ is the Immirzi parameter coming from the simplicity constraints, m, n are magnetic indices m, n = −j, −j + 1 . . . , j − 1, j, and g w ∈ SL(2, C) is a group element associated to the wedge. This restriction results from the weak quantum implementation of the simplicity constraints in the EPRL spin foam model. The Y γ map is responsible of this implementation and embeds the spin j SU (2) representation in SL(2, C) as in (2). The group element g w represents the holonomy responsible for the parallel transport along the wedge from the reference system of the source tetrahedron to the reference system of the target tetrahedron. We conventionally associate the row of the representation matrix, the couple (j, m) in (2), to the target and the column, the couple (j, n) in (2), to the source. In this way, the SL(2, C) γ-simple representation matrices inherit the orientation of the 2-complex. Instead of a group element for each wedge, we prefer to use a group element for each half-edge. We replace g w → g −1 t g s where s and t are the source and target half-edges. This choice of fundamental variables guarantees that the parallel transport on a closed path in a vertex is trivial. In other words, the product of all the holonomies on the same closed path is the identity 2 , or the holonomy is flat within a single vertex. 2 Explicitly, if w 1 , w 2 and w 3 are three wedges of the same vertex we have where we have assumed that the wedges are oriented such that the target of w 1 is the source of w 2 and so on. If the orientation of one of the wedges w is the opposite we replace gw with its inverse.
We set the spin j on each edge to be the same and contract the magnetic indices m, n. At the end of this procedure, the only non-contracted magnetic indices are on boundary half-edges. We prescribe them as part of the boundary data. A common choice to describe boundary data is to contract these magnetic indices with intertwiners in the recoupling basis or with coherent intertwiners if we are interested in representing some semi-classical geometry.
We sum over all the possible spins j f associated to each closed face and we weight the contribution of the face with the dimensional factor (2j f + 1) [32] j f where the product is on all the wedges belonging to the face. On non closed-faces we assign the spin as part of boundary data. We integrate over the group element associated to each half edge using the Haar measure of SL(2, C). For each vertex one integration is redundant and we remove it to regularize the amplitude as prescribed in [33].

Graphical notation
Writing all the constituent of an EPRL spin foam amplitude can quickly get out of hand. To help us be precise and clear, we rely on a graphical notation. We introduce the various elements as we need them. We represent a unitary irreducible representation in the principal series of SL(2, C)as an oriented line. The row labels correspond to the start of the line, and the column labels to the end of the line. We indicate the argument group element in a box and decorate the line with the needed representation labels We contract two representations summing over all the magnetic indices (both j and m are magnetic numbers from the perspective of the infinite-dimensional irreducible representations of SL(2, C)) by connecting the two lines. For example, in graphical notation, the SL(2, C) representation property reads We denote the implementation of the Y γ map (2) with a blue thick line that cuts across the representation line: If we apply the Y γ map (7) to the product g 1 g 2 and use the decomposition (6), in the graphical notation we have one blue line at both ends: The (infinite) sum over two pairs of magnetic indices is implied in graphical notation, according to equation (6) . If we contract two representation lines with a Y γ map we only sum over one pair of magnetic indices: We denote with a thicker red line the sum over the spin associated to that representation j weighted by a dimensional factor (2j + 1): The (tensor) product of two representations is represented as two lines side by side. If the group element is the same we use a single box. Similarly for the Y-map, we use a single line. When we draw a box in amplitudes, we will always imply the integration with the Haar measure over the corresponding SL(2, C) group element:

An example: writing the ∆ 4 amplitude
With the general recipe discussed in this Section and the corresponding graphical representation, we write the ∆ 4 spin foam amplitude associated with the 2-complex in Figure 2. We also inherit the naming convention from the 2-complex. We report the amplitude in Equation 12.
To assign a unique name to all the SL(2, C) group elements, we denoted as g andg the two group elements associated with the same (bulk) edge but belonging to different vertices. We used a small abuse of notation in writing (12). Some group elements appears as their inverse. To represent them as a single box we opted to not distinguish them. However, following our conventions, the group element in the matrix element of a target half-edge appears always as its inverse. For example, the half edge 1234 is the source of 234 and the target of 134. The group element g 1234 appears as D γj234,j234 (g 1234 ) and D γj134,j134 (g −1 1234 ). As mentioned above , we contracted all the boundary magnetic indices with four valent intertwiners (12 in total) as part of the prescription of the boundary data. We chose the same recoupling basis on each of them and kept the label generic for the moment (i e with e a quadruple identifying a boundary tetrahedron).
We highlighted in red the bulk face (123), dual to the triangle 123 in the ∆ 4 triangulation 1. According to equation (10) , we are implying a summation over the spin j 123 assigned to it weighted by a dimensional factor 2j 123 + 1. As part of the boundary data, we also prescribed all the spins associated with the boundary faces. We keep them generic for the moment (j f with f a triple identifying a boundary triangle).
We regularized the amplitude removing one SL(2, C) integration for each vertex as discussed above . In (12) we indicate the removed integrals with a white box. This choice is arbitrary, and the amplitude value is independent of this choice. However, we can use this arbitrariness to simplify the numerical computation (see Section 7) by making the symmetric choice. The integral removed is always opposite to the two bulk half edges and the bulk edge (123).

How-To divide the EPRL transition amplitudes
Approaching the calculation of the full amplitude is an arduous task. The group matrix elements in unitary representations are highly oscillating functions. The integrals are group integrals over many copies (sixteen in our example) of six-dimensional non-compact groups. We divide the transition amplitude into smaller and more manageable components and compute them serialized. This approach is the most advantageous if your goal is to obtain a number from the computation of a transition amplitude. However, this could be suboptimal for semiclassical calculation due to the number of components. Without any loss of generality, we insert a resolution of the identity over the intertwiner space between every two vertices of (12).g We rewrite the resolution of the identity over the intertwiner space as an integral over SU (2) of four matrix elements (56). We commute the SU (2) integral with the Y γ map and bring the SU (2) group element in the Finally, we use the invariance property of the SL(2, C) Haar measure to reabsorb the SU (2) group element with a change of variable, obtaining the original spin foam edge.

An example: decomposing the ∆ 4 amplitude
If we divide the ∆ 4 amplitude (12) inserting 4 resolutions of the identity (each one between two different vertices), the latter decomposes into a linear combination of the product of four amplitudes. That is, one per vertex. These amplitudes are commonly known as vertex amplitudes. Using the graphical representation, we write the full ∆ 4 transition amplitude as: . (16) By separating the vertices as in (16), we have transformed the problem of calculating the full amplitude into the computation of the single building blocks: the vertex amplitudes.

How-To compute the EPRL vertex amplitudes
In this section, we will focus on contributions local at the vertices. For concreteness, we model the definition of the vertex amplitude on the (12345) vertex in the example (12).
In (17) we contracted the magnetic indices of the bulk edges (1234) and (1235) with two intertwiners, labelled by i 1234 and i 1235 . We will see in the next section why this choice is natural, and we are not losing any generality. Remember that we regularized the amplitude by fixing the group element g 1245 = 1, graphically denoting such element by leaving it blank.
Consider the contribution from the wedge (234). We use the representation property to separate the matrix elements corresponding to the two group elements.
The inverse g 2345 is due to the orientation of the wedge (234) and the conventions we are adopting.
To help the reader remember about the extra summation introduced by the representation property, we wrote spin l 234 even if we are summing over it. This summation is bounded from below by j 234 but is unbounded from above. It is a consequence of the non-compactness of the group (all unitary irreducible representations are infinite-dimensional). Each group element appears as the argument of four matrix elements.
On the face (124) there is no sum over the spin l 124 , as a consequence of the regularization choice g 1245 = 1 and the presence of the Y γ map on the half-edge (1245). We parametrize each Lorentz transformation (SL(2, C) group element) with an arbitrary rotation (SU (2) group element) followed by a boost in a conventional direction (the 3 direction in our case) and another arbitrary rotation: the Cartan parametrization (68) of SL(2, C). The representation matrices decompose as (68) and the Haar measure factorizes as in (69). We divide the contribution of the integral on the half-edge (1234) to the amplitude into The matrix elements of u 1234 and v 1234 are SU (2) matrix elements elements (71). We represent the integral over the rapidity r 1234 of the product of four reduced matrix elements (73) as a white oval. We wrote the arguments explicitly to help the reader to visualize the parametrization. In the following, we will omit the name of redundant integration variables.
The Y map commutes with SU (2) group elements. Therefore, we move it next to the rapidity integral. We perform the integrals over SU (2) (58) in terms of (4jm) symbols. The contribution to the amplitude from the half-edge (20) is where the thicker red line represent a summation over the corresponding label weighted by a dimensional factor as in (10). The contraction of two (4jm) symbols obey the orthogonality condition (60) and allows us to remove the summation over i ′ where the phase (−1) 2j234 is a consequence of the different orientation of the link (46). We define the booster functions B γ 4 as the result of the integral = p1,p2,p3,p4 .
The booster functions were first introduced in [34], numerically computed in [12,35], analytically evaluated in terms of complex gamma functions [36,37], and they have an interesting geometrical interpretation in terms of boosted tetrahedra [38]. The booster functions encode how the EPRL model imposes the quantum simplicity constraints and depend on the Immirzi parameter γ. Note that, in the definition (23), we dropped the information on the orientation of the faces. The orientation of the faces in the booster function is irrelevant. The effect of orientation change of the (4jm) symbols cancels exactly the effect of orientation change of the reduced density matrices of SL(2, C), as we discuss in Appendix B. Using this definition, we can write (19) in terms of the booster functions as: .
We compute the contribution to the amplitude from all other half edges of (12345) with the same prescription. The (4jm) symbols in (24) contracts among themselves and form a {15j} symbol of the first kind (63).
The sum over spins l f are only bounded from below (e.g. l 123 ≥ j 123 ) and the intertwiners k e are bounded by the triangular inequalities of the (4jm) symbols. To complete the calculation we recognize the SU (2)invariant as a canonical {15j} symbol of the first kind (63). The vertex amplitude is In general one need to change the orientation of some links to obtain the canonical {15j} symbol using (46) to compute the relative phase.
We rewrote the vertex amplitude (17) as a combination of a canonical {15j} symbol weighted by four booster functions.

How-To calculate numbers
In the previous Section, we completed the formal evaluation of the amplitude. If we are satisfied with the expression (16) we can stop here. A few more steps are needed if we want to translate it into a number. We decompose each vertex amplitude in (16) as in (25). By doing so, we finally write the ∆ 4 transition amplitude in the appropriate form for a numerical evaluation: In this paper, we rely on the library sl2cfoam-next to perform the numerical evaluation of the EPRL spin foam amplitude. The code discussed in this Section is available in the repository in the form of notebooks [39].

Historical overview
The development of a library for the numerical computation of the Lorentzian EPRL 4-simplex vertex amplitude started with sl2cfoam [40]. The library is coded in C and is based on the decomposition of the vertex amplitude (26) in terms of booster functions. We refer to the original paper [12] for a detailed discussion of the library's performances, accuracy, and memory management. The library computes the SU (2) invariant symbols using wigxjpf [41]. The invariants are stored efficiently in custom hash tables based on khash [42] that takes into account their symmetry properties.
sl2cfoam computes the booster functions performing a numerical integration of the boost matrix elements (23). The integrand is rewritten as a finite sum of exponentials with complex coefficients to tame its highly oscillating behaviour. One obtains the booster function from the interference of many exponential integrals done with the trapezoidal rule. In order to reach enough numerical precision, the authors employed arbitrary precision floating point computations with the GNU libraries GMP [43], MPFR [44] and MPC [45].
The library was used to explore the numerical properties of the EPRL vertex [3,11,46,14]. The need for a much more efficient and accurate code immediately became clear, as the computational time for more complex amplitudes was definitely out of reach.
Recently sl2cfoam-next [47], the evolution of sl2cfoam, has been released. Also the new library is written in C, but it has an optional julia interface [48] which hugely simplifies its usage. Although sl2cfoam-next computes the Lorentzian EPRL vertex amplitude in the form (26), it introduces several ideas and techniques borrowed from High-Performance Computing and tensor networks. Therefore, with respect to the original version, it represents a significant improvement in performance and precision. We refer to [35] for its complete description and several usage examples.
The numerical integration of the booster functions is performed with the Gauss-Kronrod quadrature method after a weighted sub-intervals decomposition of the integration range. Also, for technical reasons, the γ-simple unitary irreducible representations in the principal series of SL(2, C) are slightly different from (2) as it uses D γ(j+1)j instead of D γjj .
The huge number of sums and products involved in the expression (17) is performed with optimized routines for multidimensional arrays multiplications (we will refer to them loosely as tensors in the informatics sense), such as BLAS [49] and MKL. For the description of the CPU parallelization scheme adopted, we refer to [35]. It has been recently introduced the possibility to offload tensor contractions to the GPU and parallelize them over the GPU cores with the CUDA platform [50] by using the julia package CUDA.jl [51,52].

Introducing the cut-off
The vertex amplitude (26) is made of three distinct elements: the {15j} symbol, the booster functions and the combination of two together. The formula (26) is exact and sl2cfoam-next can compute its constituents to very high numerical precision. However, the sums over the spins l f , that appear in (26) due to the split of the representation matrix elements on the wedges, are bounded from below but not from above. This means that in order to extract a number from (26) we need to make an approximation and cut-off the 6 unbounded sums in the vertex amplitude. While unbounded the sums are convergent because the vertex amplitude is finite [33]. Therefore, in principle, it is possible to find a cut-off large enough to capture the value of the amplitude with the desired precision.
The library sl2cfoam-next implement an homogeneous cut-off ∆l on all the unbounded summations. We replace the sums Unfortunately, we do not have a prescription to find the optimal value of ∆l. Numerical explorations show that it depends on the details of boundary data, such as the face spins j f and the Barbero-Immirzi parameter γ. At the moment, the best consolidated strategy is to set ∆l as large as possible and estimate the error by studying the value of the amplitude A ∆4 (∆l) as a function of the cut-off.
Recently [13] introduced an extrapolation scheme to overcome the enormous computational cost represented by indefinitely increasing ∆l. The extrapolation algorithm was used to calculate the self-energy spinfoam amplitude (see [53,46]), which is a divergent amplitude since it contains a bubble. Furthermore, we mention that the implementation of Markov Chain Monte Carlo methods in the study of spinfoams based on the techniques discussed in this paper is in progress [25].

Using the sl2cfoam-next
Before any calculation we need to import the sl2cfoam-next library and initialize it. In the blocks of code of this Section, we will imply that the library is correctly initialized first, and we omit the following code.

✝ ✆
We set the value of the Immirzi parameter to 1.2 (for historical reasons, any value is equally possible). We define a data folder that is used both to look for the fastwigxj tables and to store the computed data optionally. In this way, we avoid recomputing the same vertex amplitude for a second time. We refer to the documentation of sl2cfoam-next and the accompanying paper [35] for a detailed description of all the setup options.

Computing one vertex
We find very valuable to dedicate this paragraph to show how to use the julia frontend of sl2cfoam-next to compute the EPRL vertex amplitude. We use the amplitude (26) as reference. We provide some jupyter notebooks 3 in the Git repository [39], for interactive usage examples that the reader can compile and execute.
In the Listing 2 we show an essential schematic representation of the code.

✝ ✆
We are omitting the initialization code in Listing 1. In lines 1-2, we specify the boundary data (all the spins equal to 1) and the cut-off ∆l = 15. In line 3, we compute the amplitude. The function vertex_compute returns a tensor with five indices, one per intertwiner, computing the vertex amplitude (25) (without any phase) for all possible values of boundary intertwiners. In [39] we show how to compute a restricted range of boundary intertwiners. The @time macro is used for logging purposes, tracking the computational time and memory usage. At fixed boundary spins and Immirzi parameter, the computation time depends on several parameters such as the value of the cut-off ∆l and the accuracy level at which the library is set. With the parameters specified in Listing 1, the first time that line 3 of Listing 2 is run takes 158 seconds. The computation time decreases exponentially by selecting a lower cut-off ∆l. We tested this code on a laptop with Intel(R) Core(TM) i7-10750H 2.60GHz processor. The library distributes the workload on the available cores, according to the parallelization scheme discussed in [35]. If we store the required booster functions during the first computation, the second time we run the script takes 3.2 seconds. It is the time to compute, sum, and contract all the required {15j} symbols in the expression (25). Finally, If we store the full vertex amplitude, the computation time is negligible since nothing is calculated, and we retrieve the value from memory.
If we are interested in one single vertex amplitude this is all we need to do.

An example: computing the ∆ 4 amplitude with sl2cfoam-next
We split the computation of the amplitude (27) into two steps. First, we compute and save the value of all the necessary vertices. Then, we contract the required vertices to calculate the ∆ 4 amplitude. For simplicity, we fix all boundary spins j equal to 1. The bulk spin j 123 assumes values from 0 to 3j, while bulk intertwiners i 1234 , i 1235 , i 1236 , and i 1237 assume values compatible with triangular inequalities. With the regularization choices we did, the vertex amplitudes are fully symmetric. That is, the bulk spin and bulk intertwiners always appear in the same position in each of the four vertices. Therefore, it is sufficient to compute only a single vertex amplitude for all the possible values of spins and intertwiners. To keep the computational time reasonable, we fix the cut-off ∆l to 15. We analyze the dependence of the amplitude on this cut-off in the next step. @save "$(vertex_path)/j_bulk_$(j_bulk)_fulltensor.jld2" vertex 13 end

✝ ✆
In lines 2-3, we set all boundary spins equal to 1 and the cut-off ∆l = 15. In lines 4-6, we create the directory path to organize the files containing the computed amplitudes. In line 7, we define the range of the bulk spin, and from line 8, we loop over it. In lines 9-12, we assign the vertex amplitude's spins, compute the vertex amplitude, and save it for later use. Notice that we are computing the fulltensor vertex amplitude, namely for all the possible values of boundary intertwiners. This ensures that the ∆ 4 amplitude can be calculated for any combination of the latter. Finally we compute the whole amplitude (27) by assembling all the vertices. One of the main advantages of collecting the vertex amplitudes in multidimensional arrays is that there are very efficient methods to multiply (or "contract") the latter. For the application we discuss in this work it is unnecessary to improve upon a for loop, but julia offers the possibility to perform contractions in a wonderfully efficient and simple way, possibly using the GPU. See for example the method contract, provided in sl2cfoam-next to contract vertices with coherent boundary states. Alternatively, there are packages such as LoopVectorization.jl (see [54] for an example) or libraries like ITensor [55]. In the following block of code, we are assuming that all the variables defined in Listing 3 are available.

✝ ✆
In lines 1-2, we define the boundary intertwiners. For simplicity, we pick them all equal to 2, but any other choice is also possible. Notice that in julia the vector's index starts from 1. Therefore, we shift its value. In line 3, we initialize the variable that will contain the amplitude. From line 4, we loop over all the possible values of the bulk spin. In lines 5-6, we load the precomputed amplitude. In line 7, we initialize the variable to store the partial amplitude. The partial amplitude is the quantity in (27) at fixed value of the bulk spin j 123 . From lines 8 to 12, we sum over the bulk intertwiners the product of the four vertex amplitudes. In line 13, we add the dimensional factor to the full amplitude value, that we display in line 16.

Results and extrapolation
We summarize the result of our calculation in Table 1 and Figure 3.  Table 1: Numerical values of the amplitude A∆ 4 (∆l) in function of the cut-off. Looking at the plot in Figure 3 seems reasonable to deduce that by increasing the cut-off ∆l, the value of the amplitude grows and (asymptotically) converges to the value of the amplitude. In the first numerical works based on sl2cfoam and st2cfoam-next [11,3] the amplitude was approximated using the value with largest available cut-off. However, we have way more information (convergence, trends, speed). Is it possible to better estimate the amplitude with what we have? We use series acceleration techniques. We reorganize the sums in the amplitude such that it takes the form where a 0 is the amplitude with vanishing cut-off ∆l = 0 (also called simplified model in [34]), a 1 encodes all the terms in the various sums of A ∆4 that appears in the amplitude cut-off ∆l = 1 but are not in a 0 , and so on. The whole amplitude is recovered in the limit for infinite cutoff of (29). While recast in this form, we can apply techniques to estimate the value of numerical convergent series like the one in Appendix C. A similar approach was attempted in [13], and here we improve it and clarify it. The technique is analog to the Aitken delta-squared process [66] applied the the succession of the partial sum (29).
Since the amplitude is finite, the infinite cutoff limit of (29) exists, and the series defined in this way is convergent. We will assume that the ratios a n /a n−1 are increasing (from a certain point onward). This assumption is backed up by numerical evidence (up to the available cutoff). The lower bound estimate in (83) specialized for the series (29) is where we specified the largest maximum value of the cut-off we computed, which is ∆l = 15. The estimate (30) is significantly different from A ∆4 (15) and does not require any additional calculation or resources. The lower bound (30) is analogous to the approximation we can obtain with the Aitken's delta-squared process. With (83) we also obtain an upper bound to the amplitude.
where L = lim ∆l→∞ a n /a n−1 which we estimate numerically with a inverse power law fit as in the example in Appendix C. We stress that the validity of this upper bound needs to be taken with a grain of salt since approximating the value of L can falsify the inequality in (31). Summarizing, We can plot the amplitude together with the bound obtained from (30) and (31) to appreciate the improvement to the rough estimate A ∆4 (15). Computational resources are precious. Up to this point, all the calculations we proposed can be done on a standard laptop. How much can we improve the estimate by increasing the cut-off using High Performance Computing? We used Compute Canada's Narval cluster to increase the cut-off from 15 to 25. The computation was distributed on 80 tasks with 10 CPUs per task, requiring about 8 minutes. The script used can be found in [39]. These were the resources we could employ in this project. There is still a big room for easy improvement.  Table 2: Numerical values of the amplitude A∆ 4 (∆l) in function of the cut-off.
Repeating the estimate process, we find new upper and lower bounds.
The lower bound is marginally improved, as expected by comparing the numerical values in Table 2 to the ones in Table 1. However, the improvement on the upper bound is significant. Having more points to extrapolate the limit of the ratios L is essential. The calculation we proposed is not limited to our choice of boundary data. Using the same technique and adapting the code we compute also the value of the A ∆4 amplitude for boundary intertwiners i b = 1, 0 and for other values of the Immirzi parameter γ = 1, 0.1. We summarize the results in the Table 3 A ∆4 (25)  Table 3: Numerical calculation of the amplitude A∆ 4 with different boundary data. The boundary spins are all equal j = 1 and the cut-off ∆l = 25.

Acknowledgments
This work was made possible through the support of the FQXi Grant FQXi-RFP-1818 and of the ID# 61466 grant from the John Templeton Foundation, as part of the "The Quantum Information Structure of Spacetime (QISS)" Project (qiss.fr). This work was also supported by the Natural Science and Engineering Council of Canada (NSERC) through the Discovery Grant "Loop Quantum Gravity: from Computation to Phenomenology". We also acknowledge the Shared Hierarchical Academic Research Computing Network (SHARCNET) and Compute Canada (www.computecanada.ca) for granting access to their high-performance computing resources.
We thank F. Gozzini for very insightful comments on the numeric section of our first draft. We acknowledge the Anishinaabek, Haudenosaunee, Lūnaapéewak and Attawandaron peoples, on whose traditional lands Western University is located.
The group is homomorphic to the rotation group SO(3) and is generated by the angular momentum algebra L i with i = 1, 2, 3 satisfying the commutation relations In the fundamental representation L i = σ i /2 where σ i are the standard Pauli matrices. The Casimir operator is L 2 = L · L and the unitary irreducible representations are labeled by a spin j ∈ N/2 a half-integer and are 2j + 1 dimensional. The canonical basis for these representations diagonalizes the operator L 3 In this basis the matrix elements of the group are given by the Wigner matrices D j mn (u) ≡ j, m| u |j, n .
Their explicit expression and properties can be found in [56] and we will not report them.
In this work, we compute integrals of products of SU (2) representation matrices in terms of SU (2) invariants. We will introduce the minimal amount of tools needed and the graphical method to perform the calculations. We do not want to provide a complete introduction to recoupling theory and its graphical method that are worth books and reviews on their own [57,31,58]. We use a graphical notation that is completely analogous to the one introduced in Section 4.
We associate an oriented line to each SU (2) representation matrix. We decorate the line with a spin label and a box containing the group element We contract two representations summing over the magnetic indices by connecting the two lines. We compute the integral over SU (2) using the unique invariant measure over the group (the Haar measure du [31]). The explicit form of the measure depends on the parametrization used for the group. We collect the boxes corresponding to the same group elements. In the following, we will always imply the integration over all the group elements in the boxes.
The integral of the product of two representation matrices is given by where we defined the tensor ǫ j1 m1m2 ≡ (−1) j1−m1 δ −m1m2 , the unique invariant tensor in the product of two j 1 representations. The ǫ tensor squares to and has the symmetry property ǫ j1 m1m2 = (−1) 2j1 ǫ j1 m2m1 . We use the graphical representation to write (39) as u j 1 The invariance property of the tensor ǫ j1 m1m2 means n1,n2 From (42) we can derive the property of Wigner matrices Using this property we can also perform integrals where an inverse group element appears For simplicity, we will merge the ǫ tensors with the box in the following. At first sight, it could appear as an ambiguity since one line will have a group element u in the box, while the line with the opposite orientation u −1 and there is no indication of which is which. However, we are integrating over u. Therefore, the name we give the group element is irrelevant. The important information is contained in the relative polarity: one group element is the inverse of the other. Using this convention and the square property (40), in any closed diagram, inverting the orientation of a line (without group elements) results into a phase (−1) 2j .
The integral of the product of three representation matrices is given by The tensors appearing in (47) are the Wigner (3jm) symbols, the unique invariant tensor (or three valent intertwiner ) in the tensor product of three SU (2) representations.
The (3jm) has the following symmetry properties (see [31,58,56] for an exhaustive list) and vanishes unless the selection rules are satisfied In the graphical representation (47) is where for the (3jm) symbol we read the spins in clockwise order if all the arrows are outgoing and in anti-clockwise order if all the arrows are ingoing. In the standard SU (2) graphical calculus, this is usually indicated with a sign next to the node [58,56,57]. For our calculations, this is unnecessary, and we avoid adding this extra layer of complexity. Similarly to (44) we have The (3jm) satisfy the orthogonality relation m2,m3 and are normalized to 1 m1,m2,m3 We also compute the integral of the product of four representation matrices with the tools we provided using a trick where we doubled the number of integrals inserting a delta function. The delta function can be expanded as a sum of Wigner functions [20] δ Using the graphical representation and the properties (48) and (53) we compute (56) u j 1 where we used a red thick line to imply a summation weighted by the dimensional factor (2i + 1) and we define the (4jm) symbols (invariant tensor or four valent intertwiner) as The way we grouped representations together in (56) is completely arbitrary. The definition (59) corresponds to the choice of coupling (also called recoupling basis) of the representation of spins j 1 and j 2 (ore equivalently spins j 3 and j 4 ). The orthogonality condition (54) of the (3jm) symbols imply the normalization of the (4jm) The contraction of two (4jm) symbols in different recoupling basis forms another notable SU (2) invariant called the {6j} symbol.
The {6j} symbol in terms of (3jm) symbols can be written in a canonical from as For a numerical evaluation, it is not convenient to write the {6j} symbol as in (62). It is much more efficient to rely on libraries that compute and store Wigner {6j} symbols optimally using recursion and symmetry properties, such as wigxjpf and fastwixj [41,59].
Another higher-order invariant that appears in our calculations is the irreducible {15j} symbol of the first kind (following the classification of [58]). We can write it both graphically and in terms of {6j} symbols as: It must be emphasized that the {15j} symbol (63) is not the most convenient choice from a numerical point of view. In fact, it is possible to choose the recoupling scheme in order to obtain reducible {15j} symbols (see [60] for an example), whose evaluation is much faster. However, since this aspect is not the most critical part of the performance, we prefer to have a pleasantly symmetrical symbol and sacrifice some efficiency. This also simplifies computations of spin foams transition amplitudes with many vertices, since the basis choice in the recoupling on one edge affects both vertices it connects. Therefore, choosing a symmetric {15j} symbol as in (63), we are sure that the recoupling is consistent in every vertex.
B SL(2, C) toolbox The group SL(2, C) is the group of 2×2 complex matrices with unit determinant. The group is homomorphic to the proper Lorentz group (the Lorentz group part that preserve the sign of the time component) [61,62]. The algebra of SL(2, C) is generated by spatial rotations and boosts L i and K i satisfying the commutation relations In the spinorial representation L i = σ i /2 and K i = iσ i /2 where σ i are the standard Pauli matrices. The two Casimir operators are K 2 − L 2 and K · L. The unitary irreducible representations in the principal series are labeled by ρ a real number and k a half-integer. In these representations the Casimirs assume the values The generic unitary representation (ρ, k) is infinite dimensional since the group is non-compact. However, we can decompose the representation (ρ, k) in an infinite number of SU (2) representations that diagonalize The definition of the EPRL model is based on the canonical basis of (ρ, k). In this basis we diagonalize L 2 and L 3 L 2 |ρ, k; j, m = j(j + 1) |ρ, k; j, m , L 3 |ρ, k; j, m = m |ρ, k; j, m .
with j ≥ k and m = −j, . . . , j . The Cartan parametrization [62,34] of the group SL(2, C) is given by the map where u, v ∈ SU (2), r ∈ [0, ∞) is the rapidity and σ 3 is the diagonal Pauli matrix the generator of boosts along the z axis. The Haar measure with respect to this parametrization is [62,34] dg = 1 4π sinh 2 r dr du dv .
Moreover, e r 2 σ3 is diagonal, therefore D ρ,k jmln (e rσ3 ) = δ aa ′ D ρ,k jala (e rσ3 ) ≡ δ aa ′ d ρ,k jla (r) where d ρ,k jla are called reduced matrix elements of SL(2, C). Summarizing The expression for d ρ,k jlm (r) was given in [63,64,65,62,34] where 2 F 1 is the Gauss hypergeometric function. The phase used in (73) is the same introduced in [34], which ensures the reality of the booster function (23). The reduced matrix elements (73) satisfy the following relation: As a consequence, the matrices (71) have the property: We can write the SL(2, C) matrix elements of g −1 as: where in the first equality we used (75) (in addition to the SL(2, C) irrep properties) and in the second one (72). Since there are no phases depending on the summed index, we conclude that the orientation of the (4jm) spins in the booster function (23) is irrelevant. This justifies the fact that we draw the latter without arrows.

C Approximation of a convergent series
In this appendix, we provide further details on the extrapolation scheme used in Section 7. This is analogous to the more general Aitken's delta-squared process [66], which accelerate the rate of convergence of a sequence providing a good approximation technique . Consider the series S = ∞ n a n and cut-offed sum S N = N n a n . By definition the series is the limit of S N for infinite cut-off S = lim N →∞ S N = lim N →∞ N n a n = ∞ n a n .
Suppose that the sequence a n is positive and, from a certain point onwards, increasing such that where the ratios increase to L. The series S is convergent by the ratio test since the ratios are increasing: Hence we have a N +1 = a N aN+1 aN > a N c N , a N +2 > a N +1 c N > a N c 2 N , and in general a N +m > a N c m N for m > 0. We can provide a bound on the series observing that Similarly, we have that a k a k−1 < L ∀k > N by definition of L and monotonicity of the ratios.
Summarizing, we have an estimate from above and below of the value of the series as If the ratios are decreasing instead of increasing, we obtain an estimate analog to (82) but with inequalities operators inverted. Let's focus on (82) since it is the case relevant for the cut-off approximation presented in Section 7. We rewrite (82) in terms of cut-offed sums as Often, in real-world physical applications, the analytical expression of a n is very complicated. We cannot compute S but can still calculate the cut-offed sums S N with N as large as our numerical computational resources allow. What is the best approximation of S we can find? We will assume that we know S is convergent (so that the question is well-posed) and that the ratios c N increase. We want to use the inequalities (83). We can numerically compute the left-hand side of the inequality. What about the righthand side? The convergence of S ensure that lim N →∞ c N = L exists, however in general we cannot compute the value of L. At the moment, there is no clear strategy on how to compute L. In this work, we will consider two possibilities. None of them is optimal, and we leave improvements to future work. Notice that no matter what approximation we decide to adopt to compute L the right inequality of (83) will not hold anymore.
For example, we can approximate L with the largest available ratio L ≈ c N . In particular, if we insist and substitute in (83) the approximation L ≈ c N , the right quantity becomes equal to the left one. We have to content ourselves with a lower bound estimate of the amplitude given The estimate (84) is very similar to the strategy used in [15,13]. We clarified that it is a lower bound. Another possibility is to use the sequence of ratios c N computed numerically to estimate the value of L. This extrapolation is slightly dangerous since its accuracy depends on how large we can take the cut-off N . Of course, this is in addition to the lower bound (84).
This series is exactly summable in terms of the log function. Nevertheless, we want to approximate the series pretending not to know how to sum it, ignoring the fact that any analytical calculation is straightforward, an relying only on numerical tools. Let's assume that the largest possible cut-off we have access to is N = 15.
We can compute the cut-offed sums The largest cut-offed sum is 5% off the real value while the lower bound approximation (87) is closer, being only 0.6% off. The numerical values for the ratios are We extrapolate the limit at infinity of the ratios L fitting the data using the first few terms of an inverse power law and keeping the constant term. It is a cheap and dirty way of extrapolating, and one should be more careful. However, it is more than enough for our purposes. We use Wolfram's Mathematica built-in Fit method to perform the fit and find L ≈ 0.889. If we substitute it in (83), keeping in mind the approximations we are making, we find the estimate which is 1.6% larger than the actual value. Combining the two estimates, we obtain a range for the series S ∈ [1.549, 1.583].