Efficient Tensor Network Algorithms for Spin Foam Models

The spin foam model for the SU(2) BF theory is defined by decorating a discretized manifold or its dual 2-complex $\Delta^{\!*}$ with functions constructed from SU(2) representation theory. Each face $f$ of the 2-complex is assigned a spin $j_{f}$ labeling a SU(2) unitary irreducible representation $D^{j_{f}}$, while each edge $e$ is assigned an intertwiner $i_{e}$—an invariant map between tensor products of spin representations. The combinatorial structure of a 2-complex $\Delta^{\!*}$, which is dual to a four dimensional triangulation $\Delta$ is such that each edge constitutes four faces, while each vertex is comprised of five edges and ten faces. Each vertex is dual to a 4-simplex, with its edges dual to tetrahedra and its faces are dual to triangles.

The amplitude associated with a 2-complex $\Delta^{\!*}$ is obtained by summing over all possible assignments of bulk spin representation labels and intertwiner labels of products of local amplitudes. The amplitudes generally take the form

where $A_{f}(j_{f})=2j_{f}+1$ represents the face amplitude [45] and $A_{e}(i_{e})$ representing the edge amplitude is chosen to be dimension of the intertwiner or inverse of the norm of the intertwiner [3]. In SU(2) BF model, the vertex amplitude $A_{v}$ is given by the $\{{15j}\}$-symbol (of the first kind) [46, 40], defined in terms of Wigner $\{{6j}\}$-symbols (wigner6j) as follows:

	$\displaystyle\{{15j}\}(j_{ab};{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{e}})$	$\displaystyle=$	$\displaystyle\includegraphics[valign={c},scale={0.95}]{symbol15j.pdf}=\begin{Bmatrix}\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{1}&j_{13}&\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{3}&j_{35}&\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{5}\\ j_{12}&j_{23}&j_{34}&j_{45}&j_{15}\\ j_{25}&\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{2}&j_{24}&\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{4}&j_{14}\end{Bmatrix}$		(2)
		$\displaystyle=$	$\displaystyle(-1)^{\sum_{m}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{m}}+\sum_{n<m}j_{mn}}\sum_{x}(2x+1)\,\begin{Bmatrix}\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{1}&j_{25}&x\\ \color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{2}&j_{13}&j_{12}\end{Bmatrix}\begin{Bmatrix}\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{2}&j_{13}&x\\ \color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{3}&j_{24}&j_{23}\end{Bmatrix}\begin{Bmatrix}\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{3}&j_{24}&x\\ \color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{4}&j_{35}&j_{34}\end{Bmatrix}\times\quad$
			$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\times\begin{Bmatrix}\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{4}&j_{35}&x\\ \color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{5}&j_{14}&j_{45}\end{Bmatrix}\begin{Bmatrix}\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{5}&j_{14}&x\\ \color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{1}&j_{25}&j_{15}\end{Bmatrix}$

where the summation variable $x$, referred to as the virtual spin, ranges from $\max(|i_{1}-j_{25}|,|i_{2}-j_{13}|,|i_{3}-j_{24}|,|i_{4}-j_{24}|,|i_{5}-j_{24}|)$ to $\min(i_{1}+j_{25},i_{2}+j_{13},i_{3}+j_{24},i_{4}+j_{35},i_{5}+j_{14})$. We have followed the conventions and orientations used in [47] for the definition of the $\{{15j}\}$-symbol of the first kind in Equation (2). The intertwiners¹¹1The intertwiner labels are highlighted in red color in Equation (2). assigned to the edges are denoted by $i_{a}$, and the spins $j_{ab}$ are associated with the faces dual to the triangles.

The spin foam amplitudes (1) can be represented in a coherent basis by introducing coherent states for the boundary links [48, 39]. Each intertwiner of a boundary node is represented in the coherent state basis by group averaging the tensor product of the four coherent states associated with the node. A coherent state associated with a boundary link is characterized by a spin $j_{ab}$ and a unit normal vector ${\bf n}_{ab}\in S^{2}$ on the 2-sphere. Considering a single vertex $v$ with five boundary edges, the vertex amplitude in the coherent basis is given by the integral expression

where $G_{a}$ is a SU(2) group element in the fundamental representation, ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\cal J}:\mathbb{C}^{2}\rightarrow\mathbb{C}^{2}$ is an anti-linear map defined by $(z_{0},z_{1})\mapsto(-\bar{z}_{1},\bar{z}_{0})$, inducing a real structure on $S^{2}$, and $\langle\,,\rangle$ denote the invariant inner product on the spin $j_{ab}$ representation. One of the group integrals associated with the five edges in the definition of the coherent vertex amplitude (3) is redundant and is therefore gauge fixed to identity. The sign factor $(-1)^{\chi}$ is determined by the graphical calculus and orientations relating the spin network diagram.

The coherent states can also be introduced for every bulk face to enable a full coherent representation of the spin foam amplitude. In the integral representation of the amplitude (1), each vertex amplitude is given by the expression (3) through an insertion of a resolution of identity of each bulk spin $j$ in terms of coherent states, given by

As a result, the spin foam amplitude for a generic 2-complex is expressed as products of coherent vertex amplitudes (3). This procedure, however, introduces numerous integration variables for the normal vectors associated with the bulk spins in addition to the group integrals. The high dimensionality of the integration variables involved makes explicit evaluation of the coherent spin foam amplitudes challenging. Even for the a single vertex, where all edges are boundary (hence no integration of the normal vectors), there are a total of $12$ integration variables for the four SU(2) group integrations. The high dimensionality, combined with the highly oscillatory nature of the integrand, renders explicit numerical integrations inefficient due to slow convergence. Conversely, the integral representation is suitable for a stationary phase approximation.

Typically, the group integrals (3) can be performed analytically. Through the Peter–Weyl theorem, functions of SU(2) can be decomposed into linear combinations of functions of the spin network basis labeled by spins and intertwiners [43]. This decomposition generally allows the group integrations to be performed, resulting in sums of intertwiner variables²²2The group integrations associated with the bulk edges are replaced by sums over bulk intertwiner variables in the spin network representation., as already expressed in the amplitude (1). The coherent vertex amplitude in the spin network basis, therefore, results in a sum over intertwiner labels given by

where ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}d_{i}}:=2i+1$ is the dimension factor of the intertwiner label $i$. The coherent $\{{4j}\}$-symbol³³3A graphical notation of the coherent-$\{{4j}\}$ symbol (see reference [49]) comes with an orientation for each spin. Flipping an orientation of spin introduces a phase factor for the spin. The anti-linear map $\cal J$ acting on the normal vectors in the definition of the coherent vertex amplitude (3) leads to the choice of the orientation of the $\{{15j}\}$-symbol and coherent-$\{{4j}\}$ symbols employed in Equation (5). $\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\scalebox{1.3}{\bf c}}_{i_{e}}$ associated with the boundary edge $e$ is explicitly given in terms of the Wigner $\{\rm 3j\}$-symbols (wigner3j) by

Here, $D^{j}_{j,m}({\bf n})=\langle j,m|{\bf n}|j,j\rangle$ are the coefficients of the Wigner-$D$ matrix (wignerDjm) in the highest weight spin basis for the SU(2) group element representation of the normal vector $\bf n$. The computation of the coherent symbols (6) has been optimized in our current implementation, as available on the GitHub repository [36]. This implementation shows significant improvements over previous methods employed in [19]. The optimization is achieved by caching the wignerDjm functions since they are independent of the intertwiner index $i$ in Equation (6). This approach allows for an efficient storage, and scalability of these coherent vectors.

In summary, using coherent states as boundary data gives the coherent representation of the SU(2) BF spin foam amplitude for a 2-complex expressed in the intertwiner basis as

where $n$ represent the boundary edges (nodes) and $\chi$ depends only on the boundary spin representation labels.

Numerical computations of spin foam models primarily focus on evaluating the coherent amplitudes, such as those represented by Equation (7) in the spin-network representation. These amplitudes are expressed as sums over bulk spins and intertwiner labels. To compute the coherent amplitudes, coherent boundary data are assigned to the faces associated with each boundary edge of the 2-complex. The state-of-the-art libraries sl2cfoam [5] and its improved version sl2cfoam-next [6] are designed to perform explicit computations for Lorentzian EPRL coherent amplitudes. These libraries can also be used compute SU(2) BF coherent amplitudes. Furthermore, work in [19] utilized a numerical code in the Julia programming language to compute SU(2) BF coherent vertex amplitudes. The Lorentzian EPRL and SU(2) BF models are closely related; the boundary data in both spin foam models are given by SU(2) coherent states. Additionally, the EPRL vertex amplitudes can be expressed as an infinite, but convergent, summation over the SU(2) $\{{15j}\}$-symbols for auxiliary labels.

Various well-optimized libraries across different programming languages facilitate efficient computations of SU(2) invariants, such as Clebsch-Gordan coefficients, Wigner $\{{3j}\}$-symbols, and $\{{6j}\}$-symbols. For example, sl2cfoam and sl2cfoam-next make use of WIGXJPF and FASTWIGXJ packages within the C programming language. Similarly, the Julia programming language offers efficient packages like WignerSymbols.jl for computing SU(2) invariants, including wigner3j and wigner6j functions. These packages enable efficient computation of the $\{{15j}\}$-symbol (2), expressed as a sum of products of $\{{6j}\}$-symbols. Currently, in existing numerical implementations, each vertex amplitude represented by the $\{{15j}\}$-symbol of the first kind is computed as a 5-valent tensor in its intertwiner indices, as depicted in (8). The sum over the intertwiner indices is then performed as tensor contractions between the vertex amplitudes and the coherent vectors associated with the components of a fixed 2-complex.

However, a significant drawback of computing the $\{{15j}\}$-symbol as a tensor lies in the sheer computational complexity involved. For larger representation labels, the size of the intertwiner labels grows, and with it, the size of the $\{{15j}\}$-tensor increases exponentially, leading to severe scalability issues. For instance, given $\{{15j}\}$-tensor with uniform spins $j_{ab}=j$, each intertwiner index is of dimension $d_{j}=(2j+1)$. Therefore, the $\{{15j}\}$-tensor with equal spins has a total size of $d_{j}^{5}$ elements. For a large spin $j\gg 1$, this exponential growth in tensor size highlights the memory-intensive nature for initializing and storing them.

Moreover, the process of contracting these high-valence tensors becomes increasingly computationally expensive for larger spins due to the exponential growth in the number of operations. As an example, consider contracting a $\{{15j}\}$-tensor of uniforms spins with a vectors representing an intertwiner index, where the vector is also of dimension $d_{j}$. Such a contraction requires a computation cost of order $\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}{\cal O}(d_{j}^{5})$. When many such contractions are required, the overall computation can become prohibitively time-consuming. This generally limits the practical use of direct numerical computation methods for many vertices or for large representation labels, which are often needed to explore properties of the amplitudes and the semi-classical regime of spin foam models.

These practical limitations necessitate the exploration of more efficient numerical techniques. In the next section, we shall explore new techniques inspired by tensor network methods to address these challenges.

To illustrate this tensor network method, we first consider the coherent vertex amplitude defined in Equation (5). In the expression for the $\{{15j}\}$-symbol, each $\{{6j}\}$-symbol depends on a pair of intertwiner labels. By treating the spins $j_{ab}$ and $x$ as fixed parameters, each $\{{6j}\}$-symbol can be represented as a matrix with the pair of intertwiners as its open indices. We refer to such a matrix as the $\{{6j}\}$-intertwiner matrix, where its components are defined by:

$$w^{(j_{a},j_{b},j_{c},x)}_{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{1}\,i_{2}}:=\begin{Bmatrix}\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{1}&j_{a}&x\\ \color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{2}&j_{b}&j_{c}\end{Bmatrix}\equiv{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{1}}\,\includegraphics[scale={0.3},valign={c}]{intw_6j}\,{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{2}}\,\,.\quad\quad$$

(11)

The graphical notation in Equation (11) represents the tensor notation for the $\{{6j}\}$-intertwiner matrix $[w_{i_{1}i_{2}}]$, where the two open legs signify the intertwiner indices. Given the range of values for $i_{1}$ and $i_{2}$, the $\{{6j}\}$-intertwiner matrix can be efficiently computed for fixed values of spins $j_{a},j_{b},j_{c}$, and $x$ using existing libraries. For instance, the package WignerSymbols.jl within Julia language, has an optimized function, wigner6j, to compute the components (11) of the $\{{6j}\}$-intertwiner matrix. Additionally, known symmetries of the Wigner $\{{6j}\}$-symbol [51] can be utilized to simplify the direct implementation of these matrices. As an example, $w_{i_{1}i_{2}}$ is symmetric in its indices when either $j_{a}=j_{b}$ or $j_{c}=x$.

The expression of the $\{{15j}\}$-symbol (2) involves multiplication of $\{{6j}\}$-symbols, as components of the $\{{6j}\}$-intertwiner matrices, with no summation over the intertwiner labels shared by a pair of $\{{6j}\}$-symbols. To handle such multiplications, we introduce the following notations to represent tensors whose components are expressed as products of smaller tensors without summations over their shared indices. For example, consider a 3-valent tensor whose components are defined in terms of that of two $\{{6j}\}$-intertwiner matrices by ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}T_{i_{1}i_{2}i_{3}}}:=w_{i_{1}i_{2}}w_{i_{2}i_{3}}$, with no summation over the index $i_{2}$. This tensor can be represented using the following notation:

$$T_{i_{1}i_{2}i_{3}}=(i_{1}\,\includegraphics[scale={0.2},valign={c}]{intw_6j}\,i_{2})\,\,(i_{2}\includegraphics[scale={0.2},valign={c}]{intw_6j}\,i_{3})\,\equiv\,\begin{picture}(62.0,16.0)\put(10.0,3.0){\includegraphics[scale={0.26},valign={c}]{wig3valent} } \put(0.0,3.0){{\hbox{ i_1}}} \put(30.0,17.0){{\hbox{ i_2}}} \put(62.0,3.0){{\hbox{ i_3}},} \end{picture}$$

(12)

where the open index $i_{2}$ connects the tensor notations of the two matrices, indicating it is a shared index. Also, a vector defined by $T_{i_{1}}:=w_{i_{1}i_{1}}$ (with $i_{1}$ not summed over), which is the diagonal vector of the $\{{6j}\}$-intertwiner matrix can be represented by the notation $\includegraphics[scale={0.3},valign={c}]{dgvec}\,i_{1}$, where $i_{1}$ is its open index. These notations provide a more detailed structure of the tensor, showing the individual tensors involved.

Using the notation in (12), the $\{{15j}\}$-tensor expressed in terms of the products of components of $\{{6j}\}$-matrices⁵⁵5The representation of the $\{{15j}\}$-tensor using the $\{{6j}\}$-matrices is similar to a Matrix Product States (MPS) [29] tensor network representation. can be represented as

$$\{{15j}\}^{(J_{ab})}_{i_{1}i_{2}i_{3}i_{4}i_{5}}=\sum_{x}d_{x}\,(-1)^{i_{1}+\cdots+i_{5}}\,\begin{picture}(22.0,33.0)\put(12.0,7.0){\includegraphics[scale={0.24},rotate=180,valign={c}]{intw15j}} \put(31.0,34.0){{\hbox{ \scriptstyle i_{1}}}} \put(59.0,15.0){{\hbox{ \scriptstyle i_{2}}}} \put(50.0,-15.0){{\hbox{ \scriptstyle i_{3}}}} \put(12.0,-15.0){{\hbox{ \scriptstyle i_{4}}}} \put(5.0,16.0){{\hbox{ \scriptstyle i_{5}}}} \end{picture}\quad\quad$$

Again, note that there is no contraction or sum over all the intertwiner indices involved in the above expression. This form of the vertex amplitude is still a 5-valent tensor and is therefore, not yet computationally efficient. As discussed in section III.1, computing and storing these 5-valent tensors for evaluating coherent vertex amplitudes demands significant memory allocations, particularly for large values of the spins. To address this, we present a strategy which avoids the use these 5-valent $\{{15j}\}$-tensors, and instead rely on matrix contractions.

The strategy is based on combining coherent vectors with the $\{{6j}\}$-intertwiner matrices and reorganizing the sums involved in the coherent vertex amplitude (5). Before performing either the sum over the intertwiner labels or the virtual spin $x$ in the definition of the $\{{15j}\}$-symbol, consider the following combination

	$\displaystyle f^{(J_{12},{\bf n}_{1},\,x)}_{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{1}\,i_{2}}$	$\displaystyle:=$	$\displaystyle(-1)^{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{1}}d_{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{1}}{\scalebox{1.3}{\bf c}}_{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{1}}(j,{\bf n}_{1})\,w^{(j_{a},j_{b},j_{c},x)}_{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{1}i_{2}}$
		$\displaystyle=$	$\displaystyle(-1)^{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{1}}\left(\includegraphics[scale={0.36},valign={c}]{coh4j}\,{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{1}}\right)\left({\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{1}}\,\includegraphics[scale={0.25},valign={c}]{intw_6j}\,{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{2}}\right)$
		$\displaystyle\equiv$	$\displaystyle{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{1}}\,\includegraphics[scale={0.34},valign={c}]{coh_intw6j}\,{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{2}}.$

This combination defines a matrix $[f_{i_{1}i_{2}}]$, termed the coherent $\{{6j}\}$-intertwiner matrix, with the intertwiner labels $i_{1},i_{2}$ as its open indices. The components of the coherent $\{{6j}\}$-intertwiner matrix are given by products (without summation) of the $\{{6j}\}$-matrix $w_{i_{1}i_{2}}$ together with the coherent $\{{4j}\}$-vector ${\scalebox{1.3}{\bf c}}_{i_{1}}$, a phase factor $(-1)^{i_{1}}$, and the dimension for the intertwiner $i_{1}$. All the spins $j_{ab}$ involved in the $\{{6j}\}$-symbol and the coherent $\{{4j}\}$-symbol are collected into the parameter $J_{12}$. Additionally, the normal vectors ${\bf n}_{i}$ and the virtual spin $x$ are considered as fixed parameters. The last line of Equation (IV.1) denotes the tensor notation for the coherent $\{{6j}\}$-intertwiner matrix, where the two legs representing the intertwiner labels are its indices.

The sums over the intertwiners in the coherent vertex amplitude (5) can thus be expressed as a contraction of the coherent $\{{6j}\}$-intertwiner matrices (IV.1). By first performing the sums over the intertwiner labels, the coherent vertex amplitude (5) can be re-expressed as a sum over the virtual spin $x$ of contractions of coherent $\{{6j}\}$-intertwiner matrices as follows:

	$\displaystyle A_{v}(j,{\bf n})$	$\displaystyle=$	$\displaystyle(-1)^{\chi}\sum_{x}d_{x}\sum_{i_{1},\ldots,i_{5}}f^{(J_{12},{\bf n}_{1},\,x)}_{i_{1}\,i_{2}}\,f^{(J_{23},{\bf n}_{2},\,x)}_{i_{2}\,i_{3}}\,f^{(J_{34},{\bf n}_{3},\,x)}_{i_{3}\,i_{4}}\,f^{(J_{45},{\bf n}_{4},\,x)}_{i_{4}\,i_{5}}\,f^{(J_{51},{\bf n}_{5},\,x)}_{i_{5}\,i_{1}}$
		$\displaystyle\equiv$	$\displaystyle(-1)^{\chi}\sum_{x}d_{x}\,\,\includegraphics[scale={0.3},valign={c}]{matrix_trace}=:(-1)^{\chi}\sum_{x}d_{x}\,\,{F}(j,{\bf n},x).$

Here, ${F}(j,{\bf n},x)\in\mathbb{C}$ is defined as the trace of products/contraction of the coherent $\{{6j}\}$-intertwiner matrices for a fixed set of the spins $j_{ab},x$ and normal vectors $\bf n$ associated with the components of the vertex $v$. Its tensor network notation is depicted by the diagram in Equation (IV.1).

Thus, Equation (IV.1) provides an alternate and more efficient method to evaluate the SU(2) BF coherent vertex amplitude (5) by summing sequences of the matrix trace ${F}(j,{\bf n},x)$. This approach is computationally advantageous compared to the contraction involving the $\{{15j}\}$-tensor in Equation (10). For instance, considering a vertex with uniform spins $j_{ab}=j$ for all $a,b$, where all intertwiners have dimension $d_{j}=(2j+1)$, the computational cost of evaluating a matrix trace ${F}(j,{\bf n},x)$ is of order $\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\mathcal{O}(4\,d_{j}^{3})$. This is significantly more efficient than the contraction with a $\{{15j}\}$-tensor and a coherent vector (as discussed in Section III.1). Consequently, this method not only reduces computational complexity but also enhances scalability, making it feasible to handle larger spin representation labels.

[htb!] : SU(2) BF Coherent Vertex Amplitude

Input:
  jays: a 10-tuple of spin assignments $j_{12},j_{13},\dots,j_{45}$.

noms: sets of unit normal vectors ${\bf n}_{1},{\bf n}_{2},\dots,{\bf n}_{5}$ associated with the boundary edges.

Output: value of the coherent vertex amplitude $A_{v}\in\mathbb{C}$.

Algorithm IV.1 summarizes the necessary steps to compute SU(2) BF coherent vertex amplitudes $A_{v}$ as a sum over trace of matrices. It takes as input a set of spin assignments and unit normal vectors, and it returns the coherent vertex amplitude as a complex number. This tensor network algorithm offers an efficient method, significantly reducing memory usage and computational resources. Additionally, steps 2 through 4 within the for loop in Algorithm IV.1 can be parallelized to further optimize the computational efficiency. Algorithm IV.1 has been implemented in Julia programming language and is available on the GitHub repository [36].

Here, we adapt the tensor network algorithm to the case of Lorentzian EPRL (Engle-Pereira-Rovelli-Livine) coherent vertex amplitude. The EPRL spin foam model provides a framework for implementing a quantized version of gravity, rooted in the canonical loop quantum gravity (LQG) approach. The model is defined on a 2-complex $\Delta^{\!*}$, with each face assigned a unitary irreducible representation of the Lorentz group $\rm SL(2,\mathbb{C})$. The key idea is to implement the simplicity constraints of Plebanski formulation of gravity at the quantum level using $\gamma$-simple unitary representations of $\rm SL(2,\mathbb{C})$ [52]. Here $\gamma$ is the Barbero-Immirzi parameter.

On a 2-complex $\Delta^{\!*}$ dual to a triangulated manifold, the EPRL amplitude takes a similar form as Equation (1). The transition amplitudes are defined in terms of the $\gamma$-simple representations of $\rm SL(2,\mathbb{C})$. We refer the reader to [9, 4, 5] for details on the definition of the transition amplitudes of the EPRL model. The coherent representation of EPRL amplitudes make use SU(2) coherent states, labeled by $\{j_{\ell},{\bf n}_{\ell}\}$, as boundary data. The coherent vertex amplitude is given by

The infinite summation range over the ‘internal’ spin labels $l_{f}$ is due to the non-compactness of the gauge group. The form of the vertex amplitude arises from using the Cartan decomposition [28] of a $\rm SL(2,\mathbb{C})$ group element $g$ and an integration measure given by

where $u,v$ are $\rm SU(2)$ group elements, $\sigma_{3}$ is a Pauli matrix and $r\in[0,\infty)$ is the rapidity. The booster function $B_{4}^{\gamma}$ is defined as a one dimensional integral over the rapidity $r$ (again see [4, 5, 28] for details and definition of the Booster function). $B_{4}^{\gamma}$ encodes how the quantum simplicity constraints are implemented on a 2-complex. Furthermore, the $\{{15j}\}$-symbol is explicitly given in terms of the spin and interwiner labels by

Each booster function depends on two intertwiner indices $i_{e}$ and $k_{e}$, and hence can be represented as a matrix whose components are given by

$$b^{(J_{2},\gamma)}_{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{2}\,k_{2}}:=B_{4}^{\gamma}(j_{12},j_{23},j_{24},j_{25};\,j_{12},l_{23},l_{24},l_{25};{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{2},k_{2}})\equiv{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}i_{2}}\,\includegraphics[scale={0.38},valign={c}]{booster_mat}\,{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}k_{2}}$$

(17)

The rectangular shaped diagram in Equation (17) denotes the tensor notation for the booster functions or matrices. In the coherent amplitude (15), the sum over the intertwiner variables $i_{e},e\in\{2,\dots,5\}$ can easily be performed by contracting the booster functions with the coherent $\{{4j}\}$-vectors as

$$\sum_{\{i_{e}\}}d_{i_{e}}\,\scalebox{1.3}{\bf c}_{i_{e}}(j,{\bf n})\,b^{(J_{e},\gamma)}_{i_{e}\,k_{e}}=\includegraphics[scale={0.38},valign={c}]{booster_vec}\,{k_{e}}$$

(18)

We refer to this contraction as the coherent booster-vector, where $k_{e}$ denotes its index. Thus, the EPRL coherent vertex amplitude written as a contraction of the $\{{15j}\}$-symbol with the coherent $\{{4j}\}$ vector and booster coherent vector can be represented by the notation

$$A_{v}^{\rm EPRL}(j_{f},{\bf n})=\sum_{l_{f}=j_{f}}^{\infty}\,\,\includegraphics[scale={0.34},valign={c}]{Av_EPRL.pdf}.$$

(19)

Notice that the contraction of the EPRL coherent vertex amplitude involves four coherent booster-vectors and one coherent $\{{4j}\}$ vector. This is due to a gauge-fixing of one of the five non-compact $\rm SL(2,\mathbb{C})$ group integrals associated to the edges in its integral representation, resulting in a finite vertex amplitude. The EPRL vertex amplitude is therefore similar in structure to that of SU(2) BF model. However, due to the summation over the internal spins $l_{f}$, the tensor contractions has to be repeated many times for each possible configuration of the spins. The summation over of internal spins can be performed using approximate schemes such as truncated sums using ‘shells’ [53, 5] or acceleration operators [7, 14, 54] for faster convergence.

To optimize the computation of the coherent vertex amplitude in the EPRL model, we adopt a similar strategy akin to that used in SU(2) BF case. The coherent vertex amplitude can similarly be rewritten as a contraction of matrices by reorganizing the $\rm SL(2,\mathbb{C})$ invariant functions and the sums over the spins and intertwiners. To achieve this, we combine components of the $\{{6j}\}$-intertwiner matrix and the coherent booster-vector to give:

	$\displaystyle h^{(J_{23},\gamma,{\bf n}_{2},\,x)}_{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}k_{2}\,k_{3}}$	$\displaystyle:=$	$\displaystyle(-1)^{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}k_{2}}\sum_{i_{2}}d_{i_{2}}\,\scalebox{1.3}{\bf c}_{i_{2}}(j,{\bf n}_{2})\,b^{(J_{2},\gamma)}_{i_{2}\,{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}k_{2}}}\,w^{(j_{13},l_{24},l_{23},x)}_{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}k_{2}k_{3}}$
		$\displaystyle=$	$\displaystyle(-1)^{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}k_{2}}(\includegraphics[scale={0.37},valign={c}]{booster_vec}{\,\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}k_{2}})({\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}k_{2}\,\includegraphics[scale={0.25},valign={c}]{intw_6j}\,k_{3}})$
		$\displaystyle\equiv$	$\displaystyle{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}k_{2}}\,\includegraphics[scale={0.35},valign={c}]{coh_boostmat}\,{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}k_{3}}$

This resulting matrix is referred to as the coherent booster-matrix. The last line represents its tensor notation. Using these coherent booster-matrices, the sum over the intertwiner variables $i_{1},k_{e}$ in the coherent amplitude can be represented as a (trace of) matrix contractions. The coherent amplitude is therefore given by

	$\displaystyle A_{v}^{\rm EPRL}(j_{f},{\bf n})$	$\displaystyle=$	$\displaystyle(-1)^{\chi}\!\!\!\sum_{l_{ab}=j_{ab}}^{\infty}\!\!\sum_{x}d_{x}\!\!\!\sum_{i_{1},k_{2},\ldots,k_{5}}f^{(J_{12},{\bf n}_{1},\,x)}_{i_{1}\,k_{2}}\,h^{(J_{23},\gamma,{\bf n}_{2},\,x)}_{k_{2}\,k_{3}}\,h^{(J_{34},\gamma,{\bf n}_{3},\,x)}_{k_{3}\,k_{4}}\,h^{(J_{45},\gamma,{\bf n}_{4},\,x)}_{k_{4}\,k_{5}}\,h^{(J_{45},\gamma,{\bf n}_{3},\,x)}_{k_{5}\,i_{1}}$		(21)
		$\displaystyle=$	$\displaystyle(-1)^{\chi}\sum_{l_{ab}=j_{ab}}^{\infty}\sum_{x}d_{x}\,\includegraphics[scale={0.34},valign={c}]{Avc_eprl4}$

This formulation of the coherent vertex amplitude for the EPRL spin foam model (21) offers a more efficient alternative, and enables a faster and more scalable computations.

[htb!] EPRL Coherent Vertex Amplitude

Input:
  jays: a 10-tuple of spin assignments $j_{12},j_{13},\dots,j_{45}$

noms: sets of normal vectors ${\bf n}_{1},{\bf n}_{2},\dots,{\bf n}_{5}$ associated with the boundary edges

Output: Value of the coherent vertex amplitude $A_{v}^{\rm EPRL}\in\mathbb{C}$.

Algorithm IV.2 gives a summary of the procedure for computing the EPRL coherent vertex amplitude as a sum of sequences of matrix contractions. This algorithm leverages the matrix formulation and optimized computation techniques to efficiently evaluate the amplitude for a given set of boundary spins and normal vectors. The explicit numerical implementation of Algorithm IV.2 in Julia language is left to future work. This tensor network algorithm can also be applied within the sl2cfoam-next library where the computation of the booster functions have already been implemented.

The coherent amplitude $A_{v}$ (10) is associated with a 2-complex $v$, characterized by a single vertex with five boundary edges which is dual to a 4-simplex triangulation. To study dynamics of quantum geometries, it is essential to analyze spin foam amplitudes on a 2-complex with more than one vertex. In a generic 2-complex, with multiple vertices (see examples in Figure 1), a pair of neighbouring vertices are connected by one or more (bulk) edges. In the dual triangulation of a generic 2-complex, any 4-simplex that contains a boundary tetrahedron is dual to a vertex exhibiting a combination of bulk and boundary edges. Such a vertex is referred to as a boundary vertex, for conciseness. A bulk vertex, on the other hand, refers to a vertex with each of its five edges dual to bulk tetrahedron.

In the coherent amplitude (7), coherent data $\{j_{\ell},{\bf n}_{\ell}\}$ are assigned to the faces of boundary edges, while bulk edges are assigned spin and intertwiner data $\{j_{f},i_{e}\}$ in the spin-network basis. A boundary vertex endowed with both coherent data for the boundary edges and intertwiner data for the bulk edges is termed as a partial-coherent vertex. We will describe the amplitudes associated with these partial coherent vertices, focusing on the SU(2) BF spin foam model. Figure 2 illustrates the tensor network notations for all partial coherent vertex amplitudes, each associated with a boundary vertex.

Each of the four partial-coherent vertices in Figure 2 can be computed as a contraction of the $\{{15j}\}$-tensor with a number of coherent-$\{{4j}\}$ vectors (9), resulting in a tensor with its open legs indexed by the intertwiner labels for the the bulk edge(s) attached to the vertex. To avoid creating the $\{{15j}\}$-tensor, we shall make use of components of the $\{{6j}\}$-intertwiner matrices and the coherent $\{{6j}\}$-intertwiner matrices to rewrite these partial coherent vertices. The sums over the virtual spin and the intertwiner labels involved in the partial coherent vertices can also be reorganized such that the resulting tensors are expressed as contractions between the smaller tensors using the matrices defined in Equations (11) and (IV.1). We shall also make use of the tensor notations described in Equation (12) for the expressions.

Consider the first partial-coherent vertex with one bulk edge and four boundary edges associated with coherent data. It has one open index corresponding to the bulk intertwiner label, and hence represents a vector. Using the expressions and notations for the matrices in (11) and (IV.1), the 1-valent partial-coherent vertex amplitude can be expressed as

	$\displaystyle\includegraphics[scale={0.32},valign={c}]{pat_vertex41}\,\,i_{1}$	$\displaystyle=$	$\displaystyle(-1)^{\chi}\sum_{i_{2},\cdots,i_{5}}\,\,\{{15j}\}(j_{12},\ldots,j_{45};i_{1},\ldots,i_{5})\prod_{k=2}^{5}d_{i_{k}}{\scalebox{1.3}{\bf c}}_{i_{k}}(j,{\bf n}_{k})$		(22)
		$\displaystyle=$	$\displaystyle(-1)^{\chi}\sum_{x}d_{x}\sum_{i_{2},\ldots,i_{5}}(-1)^{i_{1}}\,w^{(J_{12},\,x)}_{i_{1}\,i_{2}}\,f^{(J_{23},{\bf n}_{2},\,x)}_{i_{2}\,i_{3}}\,f^{(J_{34},{\bf n}_{3},\,x)}_{i_{3}\,i_{4}}\,f^{(J_{45},{\bf n}_{4},\,x)}_{i_{4}\,i_{5}}\,f^{(J_{51},{\bf n}_{5},\,x)}_{i_{5}\,i_{1}}\quad$
		$\displaystyle=$	$\displaystyle(-1)^{\chi}\sum_{x}d_{x}\,(-1)^{i_{1}}\,\includegraphics[scale={0.28},valign={c},rotate=10]{cohpat_vertex41}i_{1}\,\,.$

In the second line of Equation (22), the four coherent $\{{4j}\}$-vectors associated with the boundary edges are combined with the $\{{6j}\}$-intertwiner matrices into coherent $\{{6j}\}$-intertwiner matrices $f_{ii^{\prime}}$. Just as in the case of the coherent vertex amplitude (IV.1), the sum over the intertwiner labels of the boundary edges is first performed, followed by the summation over the virtual spin $x$. Since the intertwiner label $i_{1}$ associated with the bulk edge is not summed over, the contraction results in a vector indexed by $i_{1}$. The diagram in the last line of Equation (22) represents the tensor notation of the resulting vector for fixed spins $j_{ab}$ and $x$.

The second partial-coherent vertex amplitude as shown in Figure 2 involves two bulk edges and three boundary edges with coherent data. It, therefore, represents a matrix with its indices corresponding to the intertwiner labels of the two bulk edges. The initial expression of this vertex amplitude involves a contraction of a $\{{15j}\}$-tensor with three coherent $\{{4j}\}$-vectors associated with the boundary edges. This can also be rewritten as contraction of the $\{{6j}\}$-intertwiner matrices and coherent $\{{6j}\}$-intertwiner matrices after rearranging the sums over the intertwiner labels and virtual spin as follows:

		$\displaystyle=$	$\displaystyle(-1)^{\chi}\sum_{i_{3},i_{4},i_{5}}\,\,\{{15j}\}(j_{12},\ldots,j_{45};i_{1},\ldots,i_{5})\prod_{k=3}^{5}d_{i_{k}}{\scalebox{1.3}{\bf c}}_{i_{k}}(j,{\bf n}_{k})$		(23)
		$\displaystyle=$	$\displaystyle(-1)^{\chi}\sum_{x}d_{x}\,(-1)^{i_{1}+i_{2}}\,w^{(J_{12},\,x)}_{i_{1}\,i_{2}}\sum_{i_{3},i_{4},i_{5}}w^{(J_{23},\,x)}_{i_{2}\,i_{3}}\,f^{(J_{34},{\bf n}_{3},\,x)}_{i_{3}\,i_{4}}\,f^{(J_{45},{\bf n}_{4},\,x)}_{i_{4}\,i_{5}}\,f^{(J_{51},{\bf n}_{5},\,x)}_{i_{5}\,i_{1}}\quad$
		$\displaystyle=$	$\displaystyle(-1)^{\chi}\sum_{x}d_{x}\,(-1)^{i_{1}+i_{2}}\,\begin{picture}(50.0,35.0)\put(-5.0,15.0){ \includegraphics[scale={0.28},valign={c},rotate=-36]{cohpat_vertex32} } \put(9.0,-25.0){$i_{1}$} \put(48.0,-25.0){$i_{2}$} \end{picture}\quad\quad$

The diagram in the last line of Equation (23) represent the product of the components of the matrix $w_{i_{1}i_{2}}$ with the contracted matrices involved in the sums over the boundary intertwiner labels $i_{3},i_{4},i_{5}$. These partial-coherent vertex amplitudes are required, for example, in computing the amplitude for a 2-complex with one bulk face dual to a $\Delta_{3}$ triangulation (see Figure 1). Also, a generalization to a 2-complex with one bulk face and $n$ vertices (dual to a triangulation denoted $\Delta_{n}$, for $n$ finite) implements these 2-valent parital-coherent vertex amplitudes. Permuting the boundary and bulk intertwiner labels changes the order of the contracted matrices. However, the resulting partial-coherent amplitude can still be expressed in a similar form.

The third partial-coherent vertex amplitude features three bulk edges and two boundary edges assigned with coherent data. It is a 3-valent tensor with its indices corresponding to the intertwiner labels of the bulk edges. This vertex amplitude is initially expressed as a contraction of a $\{{15j}\}$-tensor with two coherent $\{{4j}\}$-vectors associated with the boundary edges. The coherent $\{{4j}\}$-vectors and the $\{{6j}\}$-intertwiner matrices can again be combined into coherent $\{{6j}\}$-intertwiner matrices $f_{ii^{\prime}}$. The sum over the intertwiner labels gives contraction of matrices and the sum over virtual spin result in the following expression

			(24)
	$\displaystyle\begin{picture}(230.0,23.0)\put(0.0,25.0){\includegraphics[scale={0.33},valign={c}]{pat_vertex23} } \put(1.0,0.0){$i_{1}$} \put(37.0,0.0){$i_{2}$} \put(49.0,25.0){$i_{3}$} \put(62.0,25.0){$\displaystyle=(-1)^{\chi}\sum_{x}d_{x}\,(-1)^{i_{1}+i_{2}+i_{3}}\,$ } \put(183.0,11.0){\includegraphics[scale={0.28},valign={c},rotate=33]{cohpat_vertex23} } \put(193.0,0.0){$i_{1}$} \put(232.0,0.0){$i_{2}$} \put(244.0,37.0){$i_{3}$} \end{picture}\quad\quad.$

The amplitude associated with each boundary vertex of the 2-complex dual to the $T4$ triangulation, as shown in Figure 1, is given by such a 3-valent partial-coherent vertex amplitude.

Lastly, the fourth partial-coherent vertex amplitude comprises four bulk edges and one boundary edge with coherent data. It represents a 4-valent tensor with its indices corresponding to the intertwiner labels of the bulk edges. After contracting the matrices involved in the sum over the boundary intertwiner, the components of this 4-valent vertex can be represented as

In summary, each partial-coherent vertex amplitude in Figure 2 is characterized by the number of $\{{6j}\}$-intertwiner matrices and coherent $\{{6j}\}$-intertwiner matrices, corresponding to the number of bulk and boundary edges with coherent data, respectively. By reorganizing the sums over the virtual spin and intertwiner labels, the expressions for these vertices can be efficiently computed through matrix contractions, thereby enhancing memory efficiency in the computation of coherent SU(2) BF spin foam amplitudes. This new method of evaluating the partial-coherent amplitudes avoids the need to compute the $\{{15j}\}$-symbol as a 5-valent tensor, further optimizing the computation of coherent amplitudes. Although the tensor network representation of the 3-valent and 4-valent partial-coherent vertex amplitudes theoretically scales better, their performance is comparable to the computation using the $\{{15j}\}$-tensor for small representation labels. These partial-coherent vertices, together with edge amplitudes [7] are relevant for computing coherent amplitudes for generic 2-complexes with multiple boundary vertices.

As our first example, we examine a vertex $v$ with Regge boundary data, corresponding to an equilateral 4-simplex triangulation. This configuration represents the simplest and most symmetric instance of a vertex. The equilateral vertex has boundary data characterized by equal spins $j_{ab}=j$ for all faces corresponding to the areas of the dual triangles. Thus, $j$ characterizes the boundary scale. Each tetrahedron dual to an edge is also equilateral, thus for a fixed edge $a$, the unit normal vectors are explicitly given by

These set of vectors correspond to the unit normals associated with the triangular faces of each equilateral tetrahedron. They can be rotated to form a twisted spike configuration. Additionally, this set of boundary spins and normal vectors provides a consistent length geometry for an equilateral 4-simplex, with equal edge lengths given by $\ell=j(\sqrt{4}/3)$. The external dihedral angles associated with the triangles are also all equal, given by $\theta=\arccos(-\frac{1}{4})$. Appendix A provide further details for the equilateral vertex configuration. The boundary data are thus completely determined by the single boundary spin value $j$ along with the unit normal vectors (26) associated with each tetrahedron. This highly symmetric configuration significantly simplifies the computation of its coherent amplitude.

Previous numerical studies for spin foam models [19, 5, 49] have also considered the coherent amplitude for the equilateral vertex example. We will compare our results and benchmarks to those established in the earlier works, which mostly relied on contracting the $\{{15j}\}$-tensor and coherent boundary vectors. These computations require significant memory allocations and computational time, particularly for relatively large spins $(j\geq 50)$. For instance, the computation for the equilateral vertex amplitude at spin $\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}j=50$, using the $\{{15j}\}$-tensor and floating-point double precision, requires a computer with more than $\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\rm 78.3GB$ of random-access-memory⁶⁶6 Moreover, it takes more memory usage to perform the contractions of the $\{{15j}\}$-tensor with boundary coherent vectors. (RAM). Such a computation is expected to take several days even on a high-performance computing cluster with enough memory resources. There may be extra memory costs from performing contractions. These computational limits restrict the practical exploration of spin foam coherent vertex amplitudes for large spins.

In contrast, the tensor network algorithm employed here, leverages the matrix contractions to substantially reduces both computational complexity and memory usage. For instance, it now takes approximately $t\approx 1.05$ seconds (see Table 2 for more details) to compute the equilateral vertex amplitude for uniform spin $j=50$ using the improved algorithm on a consumer laptop⁷⁷7All the computations and benchmarks in this article were performed using the Julia programming language on a laptop equipped with Apple M2 Pro Chip and 16GB of RAM.. This significant reduction in computational time underscores the efficiency gained by reorganizing the computations as contractions between smaller tensors (matrices), compared to the direct computation involving the $\{{15j}\}$-tensor. Moreover, the memory-efficient nature of the new algorithm enables the exploration of coherent amplitudes at higher spin values with significantly reduced resource requirements.

Figure 3 displays the absolute value of the coherent amplitude for the equilateral vertex rescaled by a factor $j^{6}$ for the boundary spins ranging from $j=0$ to $j=200$ in steps of half integers. The rescaling is due to the power law decay $j^{-6}$ of the vertex amplitude. The plot demonstrates the capability of the new algorithms to efficiently compute coherent amplitudes across large values of the spins, which was previously impractical due to resource constraints. As shown, the coherent amplitude for the equilateral vertex configuration oscillates as a function of the spins.

The vertex amplitude for the equilateral configuration has been compared to its asymptotic formula in Appendix B, as displayed⁸⁸8The comparison shows a good agreement between the equilateral coherent vertex amplitude and its asymptotic formula with less than $\sim 1\%$ relative error for boundary spins $j\geq 100$. See Figure 8 for more details. in Figure 8. In general, the asymptotic analysis of the coherent amplitude for a vertex with Regge boundary data yields two solutions to the critical point equations of its action [55, 38]. The frequency of the oscillations is determined by the Regge action associated with the dual 4-simplex triangulation, explaining the persistence of the oscillations for large spins observed in Figure 3.

The coherent vertex amplitude for other examples of Regge geometries, such as the isosceles 4-simplex and non-regular configurations, has been studied in [15, 49]. The coherent amplitudes for the non-equilateral examples can also be efficiently computed using Algorithm IV.1. Computational results for examples of non-equilateral configurations are not presented here but are available on the GitHub repository [36].