A Landau.jl tutorial
This notebook accompanies the paper
[LD] Sebastian Mizera and Simon Telen, Landau Discriminants, arXiv: .
It illustrates the use of the Julia package Landau.jl for computing Landau discriminants. We start by adding the directory to the path and activating the environment. The latter makes sure that Julia uses the correct package versions for all auxiliary packages.
[1]:
push!(LOAD_PATH,string(pwd()))
using Pkg; Pkg.activate()
Activating environment at `~/Documents/Projects/Landau/Project.toml`
We are now ready to import the package Landau.jl. This may take a minute.
[2]:
using Landau
┌ Info: Precompiling Landau [eda78fbb-2de4-4471-838a-7157f367d5a2]
└ @ Base loading.jl:1278
In Landau.jl, a Feynman diagram is represented by a list of edges and a list of nodes. An edge between vertices \(i\) and \(j\) is represented by a 2-tuple \([i,j]\). The list of nodes indicates to which vertices the external legs are attached. For instance, for the parachute diagram in [LD, Fig. 1], the four external legs are attached to vertex 1, 1, 2 and 3 respectively.
[3]:
edges = [[3,1],[1,2],[2,3],[2,3]];
nodes = [1,1,2,3];
We can compute the Symanzik polynomials \({\cal U}_G\) and \({\cal F}_G\) for this diagram as follows.
[4]:
F, U, α, p, mm = getF(edges, nodes);
println("F = $F")
println("--------------")
println("U = $U")
F = -(α₁*α₃*α₄ + α₃*α₂*α₄)*p₃₋₄ - (α₁*α₂*α₄ + α₁*α₃*α₂ + α₁*α₃*α₄)*p₁₋₄ - (α₁*α₂*α₄ + α₁*α₃*α₂ + α₁*α₃*α₄)*p₂₋₄ - (α₁*α₂*α₄ + α₁*α₃*α₂ + α₃*α₂*α₄)*p₁₋₃ - (α₁*α₂*α₄ + α₁*α₃*α₂ + α₃*α₂*α₄)*p₂₋₃ - (α₁*mm₁ + α₂*mm₂ + α₃*mm₃ + α₄*mm₄)*(α₁*α₃ + α₁*α₄ + α₂*α₄ + α₃*α₂ + α₃*α₄)
--------------
U = α₁*α₃ + α₁*α₄ + α₂*α₄ + α₃*α₂ + α₃*α₄
The matrix p contains \(p_i \cdot p_j\) in its \((i,j)\) entry, where \(p_i\) is the momentum vector of the \(i\)-th particle. To substitute these by the masses and Mandelstam invariants, we implemented the function substitute4legs for diagrams with 4 external particles. Analogously, for diagrams with 5 dangling edges there is a function substitute5legs. These functions also has the option to use equal masses for all internal and external particles.
[5]:
F_genericmass, s, t, M, m = substitute4legs(F, p, mm);
println("F_genericmass = $(F_genericmass)")
println("-----------------------")
F_equalmass, s, t, M, m = substitute4legs(F, p, mm; equalM = true, equalm = true);
println("F_equalmass = $(F_equalmass)")
F_genericmass = (-1/2)*(-M₃ - M₄ + s)*(α₁*α₃*α₄ + α₃*α₂*α₄) + (-1/2)*(-M₂ - M₃ + t)*(α₁*α₂*α₄ + α₁*α₃*α₂ + α₃*α₂*α₄) + (-1/2)*(-M₁ - M₄ + t)*(α₁*α₂*α₄ + α₁*α₃*α₂ + α₁*α₃*α₄) + (-1/2)*(M₁ + M₃ - s - t)*(α₁*α₂*α₄ + α₁*α₃*α₂ + α₁*α₃*α₄) + (-1/2)*(M₂ + M₄ - s - t)*(α₁*α₂*α₄ + α₁*α₃*α₂ + α₃*α₂*α₄) - (m₁*α₁ + m₂*α₂ + m₃*α₃ + m₄*α₄)*(α₁*α₃ + α₁*α₄ + α₂*α₄ + α₃*α₂ + α₃*α₄)
-----------------------
F_equalmass = -(-M + (1/2)*s)*(α₁*α₃*α₄ + α₃*α₂*α₄) - (-M + (1/2)*t)*(α₁*α₂*α₄ + α₁*α₃*α₂ + α₁*α₃*α₄) - (-M + (1/2)*t)*(α₁*α₂*α₄ + α₁*α₃*α₂ + α₃*α₂*α₄) - (-M + (1/2)*(4*M - s - t))*(α₁*α₂*α₄ + α₁*α₃*α₂ + α₁*α₃*α₄) - (-M + (1/2)*(4*M - s - t))*(α₁*α₂*α₄ + α₁*α₃*α₂ + α₃*α₂*α₄) - (m*α₁ + m*α₂ + m*α₃ + m*α₄)*(α₁*α₃ + α₁*α₄ + α₂*α₄ + α₃*α₂ + α₃*α₄)
Note that the command getF also returns \({\cal U}_G\). This is computed as a byproduct since the first Symanzik polynomial appears in the definition of the second. If we want to compute just \({\cal U}_G\), we can use
[6]:
U, M, α, vtcs = getU(edges);
println("U = $U")
U = α₁*α₃ + α₁*α₄ + α₂*α₄ + α₃*α₂ + α₃*α₄
Here’s how to generate the critical point equations of \({\cal F}_G\). We load HomotopyContinuation.jl to represent these equations as a system.
[7]:
LE, α, p, mm = LandauEquations(edges, nodes);
using HomotopyContinuation
System(LE)
[7]:
System of length 4
13 variables: mm₁, mm₂, mm₃, mm₄, p₁₋₃, p₂₋₃, p₁₋₄, p₂₋₄, p₃₋₄, α₁, α₂, α₃, α₄
-mm₁*(α₁*α₃ + α₁*α₄ + α₂*α₄ + α₃*α₂ + α₃*α₄) - p₁₋₄*(α₂*α₄ + α₃*α₂ + α₃*α₄) - p₂₋₄*(α₂*α₄ + α₃*α₂ + α₃*α₄) - (α₃ + α₄)*(α₁*mm₁ + α₂*mm₂ + α₃*mm₃ + α₄*mm₄) - (α₂*α₄ + α₃*α₂)*p₁₋₃ - (α₂*α₄ + α₃*α₂)*p₂₋₃ - α₃*α₄*p₃₋₄
-mm₂*(α₁*α₃ + α₁*α₄ + α₂*α₄ + α₃*α₂ + α₃*α₄) - p₁₋₃*(α₁*α₃ + α₁*α₄ + α₃*α₄) - p₂₋₃*(α₁*α₃ + α₁*α₄ + α₃*α₄) - (α₃ + α₄)*(α₁*mm₁ + α₂*mm₂ + α₃*mm₃ + α₄*mm₄) - (α₁*α₃ + α₁*α₄)*p₁₋₄ - (α₁*α₃ + α₁*α₄)*p₂₋₄ - α₃*α₄*p₃₋₄
-mm₃*(α₁*α₃ + α₁*α₄ + α₂*α₄ + α₃*α₂ + α₃*α₄) - (α₁*α₂ + α₁*α₄)*p₁₋₄ - (α₁*α₂ + α₁*α₄)*p₂₋₄ - (α₁*α₂ + α₂*α₄)*p₁₋₃ - (α₁*α₂ + α₂*α₄)*p₂₋₃ - (α₁*α₄ + α₂*α₄)*p₃₋₄ - (α₁*mm₁ + α₂*mm₂ + α₃*mm₃ + α₄*mm₄)*(α₁ + α₂ + α₄)
-mm₄*(α₁*α₃ + α₁*α₄ + α₂*α₄ + α₃*α₂ + α₃*α₄) - (α₁*α₂ + α₁*α₃)*p₁₋₄ - (α₁*α₂ + α₁*α₃)*p₂₋₄ - (α₁*α₂ + α₃*α₂)*p₁₋₃ - (α₁*α₂ + α₃*α₂)*p₂₋₃ - (α₁*α₃ + α₃*α₂)*p₃₋₄ - (α₁*mm₁ + α₂*mm₂ + α₃*mm₃ + α₄*mm₄)*(α₁ + α₂ + α₃)
The equations we will eventually solve are an affine version of these equaitons, see [LD, Sec. 3.2]. These affine equations can be computed by the function affineLandauEquations. Notice below that they only contain 3 out of 4 Schwinger parameters (\(\alpha_4\) is set to 1) and the last equation is of the form \(yg - 1 = 0\).
[8]:
LE, y, α, p, mm = affineLandauEquations(edges,nodes);
LE, s, t, M, m = substitute4legs(LE, p, mm);
println("variables: $(variables(LE))")
println("last eq.: $(LE[end])")
variables: Variable[M₁, M₂, M₃, M₄, m₁, m₂, m₃, m₄, s, t, y, α₁, α₂, α₃]
last eq.: -1 + y*α₁*α₃*α₂*(1.0*α₁ + 1.0*α₂ + 1.0*α₃ + α₁*α₃ + α₃*α₂)
The projection of the variety defined by these 4 equations in 14 unknowns onto the space of parameters \(s,t,M_i,m_e\) is dense in the Landau discriminant. We can compute the degree of the closure of this projection as follows.
[9]:
dproj = degreeProjection(LE, [α[1:end-1];y], [s;t;M;m]);
println("degree = $dproj")
Tracking 24 paths... 100%|██████████████████████████████| Time: 0:00:09
# paths tracked: 24
# non-singular solutions (real): 6 (0)
# singular endpoints (real): 0 (0)
# total solutions (real): 6 (0)
degree = 6
The number 6 can be found in [LD, Table 1]. We can also compute samples on the Landau discriminant using the command sampleProjection, see [LD, Sec. 3.2]. The optional input npoints is a lower bound on the number of samples you would like to compute. The default is 200.
[10]:
samples, monres, AA, BB, H, singsamples = sampleProjection(LE,[α[1:end-1];y],[s;t;M;m]; npoints = 100);
println("Computed $(length(samples)) samples")
Tracking 6 paths... 100%|███████████████████████████████| Time: 0:00:04
# paths tracked: 6
# non-singular solutions (real): 6 (0)
# singular endpoints (real): 0 (0)
# total solutions (real): 6 (0)
Computed 102 samples
We now repeat the steps above, restricting to equal masses for internal and external particles. In degreeProjection and sampleProjection, we use the option findSingular = true, as it turns out that this is an example where the incidence scheme has a non-reduced component (see [LD, Remark 5]).
[11]:
LE, y, α, p, mm = affineLandauEquations(edges,nodes);
LE, s, t, M, m = substitute4legs(LE, p, mm; equalM = true, equalm = true);
dproj = degreeProjection(LE,[α[1:end-1];y], [s;t;M;m]; findSingular = true);
samples, monres, AA, BB, H, singsamples = sampleProjection(LE,[α[1:end-1];y],[s;t;M;m]; npoints = 100, findSingular = true);
println("-----------")
println("degree = $dproj")
println("Computed $(length(samples)) samples")
println("Computed $(length(singsamples)) singular samples")
Tracking 24 paths... 100%|██████████████████████████████| Time: 0:00:03
# paths tracked: 24
# non-singular solutions (real): 2 (0)
# singular endpoints (real): 6 (0)
# total solutions (real): 8 (0)
found some singular solutions
singular component is estimated to have degree >= 1
found some singular solutions
singular component is estimated to have degree >= 1
-----------
degree = 3
Computed 68 samples
Computed 34 singular samples
Both degreeProjection and sampleProjection warn that some singular solutions were found. The degree of the discriminant is 3, and from the ratio between regular and singular samples we see that it is a union of surfaces of degree 2 and 1. In particular, it is reducible. To find the equation for the degree 2 surface, we interpolate the regular samples.
[12]:
pol, N, gap = interpolate_deg(samples,2,[s;t;M;m]);
pol
[12]:
This is a polynomial with complex floating point coefficients, whose imaginary part is negligible. We convert this into a polynomial with rational coefficients via the function rat in Landau.jl.
[13]:
rat(pol)
[13]:
Similarly, for the singular samples:
[14]:
pol, N, gap = interpolate_deg(singsamples,1,[s;t;M;m]);
rat(pol)
[14]:
The product of these two polynomials (up to a scalar factor) is found when interpolating all samples simultaneously.
[15]:
pol, N, gap = interpolate_deg([samples;singsamples],3,[s;t;M;m]);
rat(pol)
[15]:
The same result is computed with the blackbox function getLandauDiscriminant, which computes the discriminants directly from the graph using the the steps illustrated above.
[16]:
Δ, par = getLandauDiscriminant(edges,nodes);
println("--------------------")
println("the discriminant is:")
Δ[1]
Tracking 24 paths... 100%|██████████████████████████████| Time: 0:00:03
# paths tracked: 24
# non-singular solutions (real): 2 (0)
# singular endpoints (real): 6 (0)
# total solutions (real): 8 (0)
found some singular solutions
singular component is estimated to have degree >= 1
found some singular solutions
singular component is estimated to have degree >= 1
found some singular solutions
singular component is estimated to have degree >= 1
Interpolating 24 samples ...
M*m^2 + (-11/19)*M^2*m + (-4/19)*m^2*s + (4/19)*M*m*s + (1/19)*M^3 + (-9/19)*m^3
--------------------
the discriminant is:
[16]:
If possible, this function returns the individual components of the discriminants. Here is the example of the double box (dbox) diagram.
[17]:
edges = [[1,2],[2,3],[3,4],[4,5],[5,6],[6,1],[3,6]]
nodes = [1,2,4,5]
Δ, par = getLandauDiscriminant(edges,nodes);
Tracking 436 paths... 100%|█████████████████████████████| Time: 0:00:04
# paths tracked: 436
# non-singular solutions (real): 8 (0)
# singular endpoints (real): 0 (0)
# total solutions (real): 8 (0)
Components have degrees [2, 4], total estimated degree is 6.
-------------------------------------------
Sampling component 1 ...
Interpolating 12 samples ...
-M*m + (1/4)*m*s + m*t + (-1/4)*s*t
-------------------------------------------
Sampling component 2 ...
Interpolating 44 samples ...
M^2*m^2 + (-2/3)*M^3*m + (1/16)*m^2*s^2 + (1/9)*m^2*t^2 + (1/144)*s^2*t^2 + (-1/2)*M*m^2*s + (-2/3)*M*m^2*t + (1/6)*M^2*m*s + (1/9)*M^2*m*t + (-1/18)*M^2*s*t + (-1/18)*m*s*t^2 + (-5/72)*m*s^2*t + (1/6)*m^2*s*t + (5/18)*M*m*s*t + (1/9)*M^4
[18]:
Δ[1]
[18]:
[19]:
Δ[2]
[19]:
Finally, we show how to compute the number of master integrals of the double box diagram. The function \(\chi\) implements exactly the code snippet shown in [LD, Sec. 5.2].
[22]:
Landau.χ(edges,nodes)
- <span class=”ansi-green-fg”>Solutions found: 227 Time: 0:00:06</span>
<span class=”ansi-blue-fg”> tracked loops (queued): 1589 (0)</span> <span class=”ansi-blue-fg”> solutions in current (last) loop: 0 (0)</span> <span class=”ansi-blue-fg”> generated loops (no change): 7 (5)</span> </pre>
- textcolor{ansi-green}{Solutions found: 227 Time: 0:00:06}
textcolor{ansi-blue}{ tracked loops (queued): 1589 (0)} textcolor{ansi-blue}{ solutions in current (last) loop: 0 (0)} textcolor{ansi-blue}{ generated loops (no change): 7 (5)} end{sphinxVerbatim}
- [32mSolutions found: 227 Time: 0:00:06[39m
[34m tracked loops (queued): 1589 (0)[39m [34m solutions in current (last) loop: 0 (0)[39m [34m generated loops (no change): 7 (5)[39m
[22]:
227
The number 227 can be found in [LD, Tab. 3].
[ ]: