In this notebook it is shown how to:
The standard reference is arXiv:2107.13952 (henceforth denoted as [F]).
An excellent additional resource to this notebook is arXiv:2202.04360 by Donà, Frisoni (henceforth denoted as [DF]), where you can find a step-by-step tutotial on how to compute spinfoam amplitudes using sl2cfoam-next from the very basics. Their companion repository walks you through the calculation of the $\Delta_4$ amplitude.
You can find the code at this Github repository. Follow the Quickstart guide for installation. Installation might be tricky, especially making the compilation works with your current instances of Intel's MKL. Therefore I recommend compiling with CBLAS/openBLAS instead.
Assuming you followed the Quickstart guide mentioned above, the following code initializes the library. Refer to [F] or the repo for the various options.
using SL2Cfoam
# Immirzi parameter, can be changed later
Immirzi = 0.2
# working folder that also contains precomputed Wigner-nj symbols
sl2c_data_folder = "/home/francesco/hd/src/data_sl2cfoam_next"
# configuration with options
sl2c_configuration = SL2Cfoam.Config(VerbosityOff, VeryHighAccuracy, 100, 0)
# call to init
SL2Cfoam.cinit(sl2c_data_folder, Immirzi, sl2c_configuration)
The library is now initialized. We can start using it.
The most basic and most important functionality of the library is to compute EPRL vertex amplitudes in the spin-intertwiner basis. The simplest input to the vertex_compute function is the list of the 10 boundary spins and an integer number of shells (homogeneous cutoff on the virtual spins from Speziale's decomposition). Refer to the source code (here in particular) for optional additional parameters.
The amplitudes are computed for all possible values of the boundary intertwiners ($\{0, 1,2\}$ in this case) and packed into a Vertex object.
# spins
js = ones(10)
# shells
Dl = 2
# compute the amplitude
v = vertex_compute(js, Dl);
@show size(v.a);
@show length(v.a);
size(v.a) = (3, 3, 3, 3, 3) length(v.a) = 243
Intertwiners can be contracted among themselves. The contract function does that (also a simple for loop would suffice) for boundary coherent states. The following code generates a coherent intertwiner (corresponding to a semiclassical tetrahedron, see [DF]) and contracts it alng the first intertwiner $i_5$.
# define the 4 spins that enter the coherent state
js_intw = ones(Spin, 4)
# compute the coherent state with random angles in (0, pi)
# angles are 4x2 matrices (theta, phi) [ 2 angles per 4 normals ]
cs5 = coherentstate_compute(js_intw, pi * rand(4,2));
# contract the vertex with the coherent state -> we get a matrix one-dimension lower
size(contract(v, cs5))
(3, 3, 3, 3)
One can also contract vertices among themselves, using e.g. a for loop (a faster option would be to use vectorization, but Julia is actually quite smart in optimizng this sort of things). An example will be provided later when computing a simple graph with internal faces.
One of the simplest calculations is to check the asymptotic scaling of the vertex amplitude. In this setting, some semiclassical coherent state is imposed on the boundary and the boundary spins are rescaled. Barret's theorem says that the amplitude will oscillate approximately as $ \exp iS_\text{Regge}$ in the limit, where $S_\text{Regge}$ is the Regge calculus action of the semiclassical geometry imposed on the boundary. See arXiv:1903.12624 for precise statements about the theoretical analysis of this setting.
Imposing random angles on the boundary tetrahedra, the amplitude is suppressed because the closure contraint is not satisfied.
The first cell calculates all the needed vertices and stores them. It takes ~3 minutes on my laptop. You can set a lower value of Dl (the number of shells, above I set Dl = 2) to speed up the calculation significantly.
# compute vertices
# all boundary spins set to λ
@time vs = [ vertex_compute(λ * ones(10), Dl) for λ in 1:10 ];
163.865426 seconds (66.50 M allocations: 2.960 GiB, 0.14% gc time, 0.12% compilation time)
Now we contract with coherent states with random angles to obtain the amplitude (a complex number for each value of $\lambda$)
ampls_rand = ComplexF64[]
# set of 5 matrices of 4 random angles (theta, phi)
# in spherical coordinates, one per each boundary tetrahedron
angles_rand = [ pi * rand(4,2) for i in 1:5 ]
for λ = 1:10
# compute 5 coherent intertwiners with random angles
css = [ coherentstate_compute(λ * js_intw, angles_rand[i]) for i in 1:5 ]
# contract the vertex with the 5 intertwiners
ampl = contract(vs[λ], css...)
# store the amplitude
push!(ampls_rand, ampl)
end
using Plots
λs = collect(1:10)
# logplot of absolute value of the coherent amplitudes
plot(λs, λs.^(-12) * abs.(ampls_rand)[1], label = "asymptotic")
scatter!(λs, abs.(ampls_rand), label = "random angles")
plot!(yaxis=:log, title="Simple asymptotic analysis", legend = :bottomleft)
We see that the amplitude is suppressed with respect to the expected power law $\lambda^{-12}$.
Imposing the angles of a geometrical 4-simplex, we can find that the amplitude is not suppressed anymore. Notice that now the calculation is almost istantaneous since the computed vertices have been stored, and now only the calculation of coherent states and contractions are left.
In the first cell we compute the angles from the normals given in arXiv:1903.12624, Table 4.
normals = [
[(0, 0, 1), (0.94, 0, -0.33), (-0.47, 0.82, -0.33), (-0.47, -0.82, -0.33)],
[(0, 0, -1), (-0.24, -0.91, 0.33), (0.91, 0.25, 0.33), (-0.67, 0.66, 0.33)],
[(-0.94, 0, 0.33), (0.24, 0.91, -0.33), (0.09, -0.66, -0.75), (0.62, -0.25, 0.75)],
[(0.47, -0.82, 0.33), (-0.91, -0.25, -0.33), (-0.09, 0.66, 0.75), (0.53, 0.41, -0.75)],
[(0.47, 0.82, 0.33), (0.67, -0.66, -0.33), (-0.62, 0.25, -0.75), (-0.53, -0.41, 0.75)]
]
angles = []
for tet in normals
angles_tet = zeros(4, 2)
for (i, n) in enumerate(tet)
theta, phi = acos(n[3]), atan(n[2], n[1])
angles_tet[i, :] = [theta, phi]
end
push!(angles, angles_tet)
end
Now we contract with coherent states with these angles to obtain the amplitude (a complex number for each value of $\lambda$)
ampls_eucl = ComplexF64[]
for λ = 1:10
css = [ coherentstate_compute(λ * js_intw, angles[t]) for t in 1:5 ]
ampl = contract(vs[λ], css...)
push!(ampls_eucl, ampl)
end
scatter!(λs, abs.(ampls_eucl), label = "4-simplex angles")
Now we see that the amplitude follows the expected power law $\lambda^{-12}$.
A finer analysis would reveal also the oscillations with the correct phase. The Lorentzian case is more computationally demanding as the boundary spins are necessarily higher in order to reconstruct a Lorentzian geometry, and a higher number of shells is needed. This is discussed thoroughly in [F].
In this section we do some simple calculations for the self-energy graph (also called bubble or melon graph in the literature), the contraction of 2 vertices along 4 intertwiners. It is the spinfoam equivalent of the simplest dominant loop correction to the propagator in QFT. I follow arXiv:2112.14781 to which I refer for many details.
The graph has 4 boundary faces/spins which I fix to $\tfrac{1}{2}$ and 6 internal faces/spins which are unconstrained. The final amplitude is the sum over all internal unbounded spins and divergent in the cutoff over the internal spins, as commonly happens in QFT when calculating loop corrections.
First, we need the function to compute a single contracted amplitude, given a choice of the 6 internal spins. For the purpose of this simple introduction we set zero shells (also called simplified model).
using LinearAlgebra
using HalfIntegers
# I set the 4 boundary spins
onehalf = half(1)
j12 = j13 = j14 = j15 = onehalf
# the following option tells vertex_compute NOT to store the amplitudes
# because a huge number of files would be created otherwise
# the boosters are stored anyway (they are much less in number)
# so that re-evaluating the function takes less time
result_return = (ret = true, store = false, store_batches = false);
function bubble(j23, j24, j25, j34, j35, j45)
# compute vertex
v = vertex_compute([j12, j13, j14, j15, j23, j24, j25, j34, j35, j45], 0; result = result_return)
# trick for contraction:
# fixing only one boundary intertwiner, Julia can contract the rest of the vertex
# by simply considering it a giant vector, thanks to the chosen memory layout
# (also for loops would work)
# the following line selects i_1 = 0 (possible values are [0, 1])
# and creates a view of the corresponding 4-dimensional submatrix
sub_array = @view v.a[:, :, :, :, 1]
# contraction is just dot product
ampl = dot(sub_array, sub_array)
# add dimensional factors and return
ampl * dim(j23) * dim(j24) * dim(j25) * dim(j34) * dim(j35) * dim(j45)
end
bubble (generic function with 1 method)
# test: notice that not all choices of spins are allowed because of triangular inequalities;
# non-allowed choices give zero
@show bubble(1/2 * ones(6)...);
@show bubble(ones(6)...);
bubble((1 / 2) * ones(6)...) = 1.3866539766959758e-8 bubble(ones(6)...) = 0.0
Now we need a function to loop over all internal spins, up to a maximum homogenous cutoff $K$. We are interested in the values at the "partial cutoffs" $k \leq K$ to study the law of divergence.
A technical digression. The library parallelizes over virtual spins of the vertex from the shells expansions using OMP. In the present calculation, it is more convenient to parallelize over the internal spins of the graph. Therefore I turn off the parallelization over shells and write a parallel Julia loop (the Julia kernel must be setup to use more than 1 thread).
# disable library automatic parallelization
SL2Cfoam.set_OMP(false)
function bubble_divergence(cutoff)
ampls = Float64[]
# this ensures concurrent access to the summed amplitude by all Julia threads
ampl = Threads.Atomic{Float64}(0.0)
# track all spins to be computed
spins_tot = 0
# loop over partial cutoffs
for k = 0:onehalf:cutoff
# generate a list of all spins to compute
spins_all_k = NTuple{6, HalfInt}[]
for j23::HalfInt = 0:onehalf:k, j24::HalfInt = 0:onehalf:k, j25::HalfInt = 0:onehalf:k,
j34::HalfInt = 0:onehalf:k, j35::HalfInt = 0:onehalf:k, j45::HalfInt = 0:onehalf:k
# skip if computed in lower partial cutoff
j23 < k && j24 < k &&
j25 < k && j34 < k &&
j35 < k && j45 < k && continue
# skip if any intertwiner range is empty
r2, _ = intertwiner_range(j12, j25, j24, j23)
r3, _ = intertwiner_range(j23, j13, j34, j35)
r4, _ = intertwiner_range(j34, j24, j14, j45)
r5, _ = intertwiner_range(j45, j35, j25, j15)
isempty(r2) && continue
isempty(r3) && continue
isempty(r4) && continue
isempty(r5) && continue
# must be computed
push!(spins_all_k, (j23, j24, j25, j34, j35, j45))
end
# save all spins to be computed
spins_tot += length(spins_all_k)
# parallel Julia loop
Threads.@threads for spins in spins_all_k
j23, j24, j25, j34, j35, j45 = spins
# calculate the amplitude
a = bubble(j23, j24, j25, j34, j35, j45)
Threads.atomic_add!(ampl, a)
end # threads
println("Amplitude at partial cutoff = $k: $(ampl[])")
flush(stdout)
push!(ampls, ampl[])
end # partial cuttofs loop
println("Done. Computed $spins_tot single bubble amplitudes.")
ampls
end
bubble_divergence (generic function with 1 method)
Now we can run a calculation. The first evaluation takes roughly a minute on my laptop to go up to maximum cutoff $K = 6$ in this simple setting. I have disabled automatic storage of the amplitude tensors as there's a huge number of them and it might stress the hard drive a bit. Neverthless the boosters are stored anyway so that subsequent evaluations will run much faster.
cutoff = 6
@time bubble_ampls = bubble_divergence(cutoff);
Amplitude at partial cutoff = 0: 0.0 Amplitude at partial cutoff = 1/2: 1.471907125568795e-7 Amplitude at partial cutoff = 1: 3.7187466189494386e-7 Amplitude at partial cutoff = 3/2: 5.66954360959087e-7 Amplitude at partial cutoff = 2: 7.203370524737356e-7 Amplitude at partial cutoff = 5/2: 8.414556417074719e-7 Amplitude at partial cutoff = 3: 9.375086939170536e-7 Amplitude at partial cutoff = 7/2: 1.0150810692032158e-6 Amplitude at partial cutoff = 4: 1.0784740900438835e-6 Amplitude at partial cutoff = 9/2: 1.1310806540703103e-6 Amplitude at partial cutoff = 5: 1.1752121000671041e-6 Amplitude at partial cutoff = 11/2: 1.2126808277755375e-6 Amplitude at partial cutoff = 6: 1.2447804956109272e-6 Done. Computed 163647 single bubble amplitudes. 67.934476 seconds (561.72 M allocations: 17.249 GiB, 0.20% gc time, 0.86% compilation time)
ks = collect(0:0.5:cutoff)
scatter(ks, bubble_ampls, legend = :none, title = "Self-energy divergence")
Of course, this simple analysis is quite inconclusive. I refer to arXiv:2112.14781 for convincing numerical evidence of a scaling that is linear in the cutoff, when taking the double limit $K, \Delta s \to \infty$. The exact scaling is not known analytically, see references therein.
Since in this example we have a summation over many variables, we can leverage Montecarlo (MC) sampling techniques. The idea is to sample over the 6-dimensional discrete space of internal spins. We use a uniform prior for this simple example. A complete calculation with many more details and a refined analysis is found in arXiv:2302.00072.
Since the lower cutoff calculations can be performed exactly in short time and the MC approximation won't work well with small sample spaces, we use it go beyond $K = 6$ calculated previously.
# function to sample the list of internal spins
# takes as input cutoff to start, cutoff to end (Kmin < k <= Kmax)
# and a maximum number of samples per partial cutoff
function sample_spins_MC(Kmin, Kmax, N)
# a dictionary indexed by the partial cutoffs
spins_sampled = Dict{HalfInt, Vector{NTuple{6, HalfInt}}}()
# loop over partial cutoffs
for k = (Kmin + onehalf):onehalf:Kmax
# generate a list of all spins to compute
spins_sampled_k = NTuple{6, HalfInt}[]
for n in 1:N
# sample six spins from the uniform distribution for this partial cutoff
j23 = half(rand(0:2*k))
j24 = half(rand(0:2*k))
j25 = half(rand(0:2*k))
j34 = half(rand(0:2*k))
j35 = half(rand(0:2*k))
j45 = half(rand(0:2*k))
# skip if computed in lower partial cutoff
j23 < k && j24 < k &&
j25 < k && j34 < k &&
j35 < k && j45 < k && continue
# skip if any intertwiner range empty
r2, _ = intertwiner_range(j12, j25, j24, j23)
r3, _ = intertwiner_range(j23, j13, j34, j35)
r4, _ = intertwiner_range(j34, j24, j14, j45)
r5, _ = intertwiner_range(j45, j35, j25, j15)
isempty(r2) && continue
isempty(r3) && continue
isempty(r4) && continue
isempty(r5) && continue
# add to sample
push!(spins_sampled_k, (j23, j24, j25, j34, j35, j45))
end
# add to dictionary
spins_sampled[k] = spins_sampled_k
end
return spins_sampled
end
sample_spins_MC (generic function with 1 method)
# calculates the divergence with MC
# takes as input the amplitude at Kmin and
# the spins sampled with the previous function
function bubble_divergence_MC(ampl_Kmin, spins_sampled, N)
ampls_MC = Float64[]
ampl_MC = ampl_Kmin
# get partial cutoffs to do
ks = sort(collect(keys(spins_sampled)))
# loop over partial cutoffs
for k in ks
# this ensures concurrent access to the array by all Julia threads
ampl_k = Threads.Atomic{Float64}(0.0)
spins_sampled_k = spins_sampled[k]
# parallel Julia loop
Threads.@threads for spins in spins_sampled_k
j23, j24, j25, j34, j35, j45 = spins
# calculate the amplitude
a = bubble(j23, j24, j25, j34, j35, j45)
Threads.atomic_add!(ampl_k, a)
end # threads
# result = MC average * volume
# the volume is the number of points in the hypercube with (2k + 1) points per edge, in 6 dimensions
# the MC average picks up the average amplitude in the subvolume contained
# in the partial cutoff layer intersected with the constraints from triangular inequalities
volume = (2k + 1)^6
ampl_MC += ampl_k[] / N * volume
println("Amplitude at partial cutoff with MC = $k: $ampl_MC")
flush(stdout)
push!(ampls_MC, ampl_MC)
end # partial cuttofs loop
ampls_MC
end
bubble_divergence_MC (generic function with 1 method)
Now I calculate the higher cutoff points with MC, using sample sizes $N = 1000, 10000$ and check the error.
# nr of samples
Nk = 1000
# sample
@time sample = sample_spins_MC(6, 8, Nk);
# compute
@time bubble_ampls_MC1 = bubble_divergence_MC(bubble_ampls[end], sample, Nk);
0.203911 seconds (312.24 k allocations: 16.133 MiB, 97.93% compilation time) Amplitude at partial cutoff with MC = 13/2: 1.3101811043349294e-6 Amplitude at partial cutoff with MC = 7: 1.3257616688436107e-6 Amplitude at partial cutoff with MC = 15/2: 1.3333444728431556e-6 Amplitude at partial cutoff with MC = 8: 1.3392526146494814e-6 2.955936 seconds (20.74 M allocations: 800.757 MiB, 0.56% gc time, 9.22% compilation time)
# nr of samples
Nk = 10000
# sample
@time sample = sample_spins_MC(6, 8, Nk);
# compute
@time bubble_ampls_MC2 = bubble_divergence_MC(bubble_ampls[end], sample, Nk);
0.050460 seconds (1.02 M allocations: 20.016 MiB) Amplitude at partial cutoff with MC = 13/2: 1.2676798331549492e-6 Amplitude at partial cutoff with MC = 7: 1.2948862739666276e-6 Amplitude at partial cutoff with MC = 15/2: 1.3241075464184536e-6 Amplitude at partial cutoff with MC = 8: 1.342256929115057e-6 29.292441 seconds (260.95 M allocations: 9.844 GiB, 0.08% gc time)
Finally I check the exact numbers running again the exact code (takes a while...)
cutoff = 8
@time bubble_ampls_extended = bubble_divergence(cutoff);
Amplitude at partial cutoff = 0: 0.0 Amplitude at partial cutoff = 1/2: 1.4719071255687948e-7 Amplitude at partial cutoff = 1: 3.718746618949438e-7 Amplitude at partial cutoff = 3/2: 5.669543609590866e-7 Amplitude at partial cutoff = 2: 7.203370524737352e-7 Amplitude at partial cutoff = 5/2: 8.414556417074715e-7 Amplitude at partial cutoff = 3: 9.375086939170531e-7 Amplitude at partial cutoff = 7/2: 1.0150810692032164e-6 Amplitude at partial cutoff = 4: 1.078474090043884e-6 Amplitude at partial cutoff = 9/2: 1.131080654070311e-6 Amplitude at partial cutoff = 5: 1.1752121000671047e-6 Amplitude at partial cutoff = 11/2: 1.2126808277755381e-6 Amplitude at partial cutoff = 6: 1.2447804956109278e-6 Amplitude at partial cutoff = 13/2: 1.2725411763330556e-6 Amplitude at partial cutoff = 7: 1.2967274438100262e-6 Amplitude at partial cutoff = 15/2: 1.317959298691506e-6 Amplitude at partial cutoff = 8: 1.3367118017663789e-6 Done. Computed 709860 single bubble amplitudes. 166.675033 seconds (1.39 G allocations: 48.236 GiB, 0.55% gc time)
Ks_higher = collect(5:0.5:cutoff)
scatter(Ks_higher, bubble_ampls_extended[11:end], label = "exact", title = "Self-energy divergence")
scatter!(Ks_higher, vcat(bubble_ampls, bubble_ampls_MC1)[11:end], label = "MC, N = 1000")
scatter!(Ks_higher, vcat(bubble_ampls, bubble_ampls_MC2)[11:end], label = "MC, N = 10000")
We see the quite good agreement of the MC calculation, which is much faster already at this low cutoff. Notice that each run of the MC sampling will produce a different plot. One can simply take the average of many runs to get a good estimate plus error.