Computerised calculation of (SuperMinus[e])(SuperPlus[e])->(SuperMinus[μ])(SuperPlus[μ])

In this Mathematica notebook, I will demonstrate how to use computers to calculate (SuperMinus[e])(SuperPlus[e])->(SuperMinus[μ])(SuperPlus[μ]). This supplements the calculation in Section 6.3 of my QFT lecture notes and uses the same notation etc. The goal of this notebook is to demonstrate that these calculations can be fully automatised, at least as long as there are no loops are involved.

This notebook runs through on my Linux laptop where I have installed all dependencies. If you want to run this on your own laptop, please install Package X by downloading this zip file and uncompressing it in your $UserBaseDirectory/Applications/ directory. Please also download the file qgraf.wl into the same directory.

To find this folder, you can execute the following code block

SystemOpen[FileNameJoin[{$UserBaseDirectory,"Applications"}]]

Installing the diagram generator is more complicated and will not be part of this example.

Load Package X and QGraf

Once you have installed Package X and qgraf.wl we can load it here. Note that this will fail if you have not installed it correct.

<<X`

Get[FileNameJoin[{$UserBaseDirectory,"Applications","qgraf.wl"}]]

Following the recipe

Step 1-4: generating the amplitude

For the first four steps

Draw all Feynman diagrams.

assign momenta to all edges and Lorentz indices to all photon vertices

pick the end of any fermion line and follow the arrow backwards, evaluating as you go

multiply in the polarisation tensors for any external photons and propagators for each internal photon

We can use the diagram generator QGraf. Note that this section of the code has been disabled to avoid accidental deletion

amplitude=RunQGraf[{"elm","elp"},{"mum","mup"},0][[;;,2]]/.line[a_,b_]:>b/._δZ|_δm->0

QGrafGetDiagrams["qgraf[{elm,elp},{mum,mup},0]"]

amplitude={(((e^2) ⟨𝓋[p2,m],(Subscript[γ, μ1]),𝓊[p1,m]⟩⊗⟨𝓊[q1,M],(Subscript[γ, μ2]),𝓋[q2,M]⟩ (Subscript[𝕘, μ1,μ2]))/(2 (m^2)+2 p1.p2))};

As you can see, the code has generated the diagram we were after using the correct spinors but has dropped a factor of i. Since we will calculate ((Abs[ℳ])^2) this would drop out anyway. Also, note that the code has picked momenta labels (Subscript[p, i]) for incoming and (Subscript[q, i]) for outgoing momenta as well as ‘random’ Lorentz indices.

Please note that Package-X uses always lower indices. This is not formally correct but simplifies the notation and it can always be mapped on the more correct notation.

Step 5-7: finishing the amplitude

Let us relabel the outgoing momenta to match the lecture notes

contract any open index

Contract[I amplitude]

if your diagram has any loops, integrate over the unconstrained momenta. If you have any internal fermion loops, calculate the trace over their gamma matrices.

Since we have no loops to integrate, we do nothing

you now have the amplitude ℳ

ℳ=Contract[I Total[amplitude]/.{q1->p3,q2->p4}]

Steps 8-10: squaring the amplitude

We can perform these operations manually or we use the function FermionSpinSum to perform steps 8-10

square the amplitude as ((Abs[ℳ])^2)=ℳ ℳ*

ℳconj=-ComplexConjugate[ℳ]

Note how we have now introduced the index (Subscript[μ, 22]) to avoid collisions!

if you calculate unpolarised scattering, sum over final-state and average over initial-state polarisations.

for fermions, this will result in a trace of γ matrices

ℳsq=FermionSpinSum[ℳ ℳconj]

We can now contract the open indices to write terms of scalar products.

ℳsqS=Simplify[Contract[ℳsq]/.𝒹->4]

We could express this in terms of s, t, u using

Simplify[ℳsqS/.MandelstamRelations[{p1,p2,p3,p4,m,m,M,M}->{s,t,u}]]

Kinematics

$Assumptions=s>4(M^2)>4(m^2)>0&&1>cos>-1&&2π>ϕ>0&&1>β>0;

If we want to calculate for example, dσ/dΩ it would be better to write this in terms of cosθ and ϕ. To do this, let us define the scalar product as

LDot[{e1_,x1_,y1_,z1_},{e2_,x2_,y2_,z2_}]:=e1 e2-x1 x2-y1 y2-z1 z2

and parameterise our momenta (here we once again use a function that is not part of the normal Mathematica suite, therefore the following block is disabled).

momenta=Thread[{p1,p2,p3,p4}->PSX2[{(ϕ/(2π)),((cos+1)/2)},s,{m,m,M,M}][[;;,{4,1,2,3}]]]

momenta={
p1->{((Sqrt[s])/2),0,0,(Sqrt[-(m^2)+(s/4)])},
p2->{((Sqrt[s])/2),0,0,-(Sqrt[-(m^2)+(s/4)])},
p3->{((Sqrt[s])/2),(((Sqrt[1-(cos^2)]) (Sqrt[s (-4 (M^2)+s)]) Cos[ϕ])/(2 (Sqrt[s]))),(((Sqrt[1-(cos^2)]) (Sqrt[s (-4 (M^2)+s)]) Sin[ϕ])/(2 (Sqrt[s]))),((cos (Sqrt[s (-4 (M^2)+s)]))/(2 (Sqrt[s])))},
p4->{((Sqrt[s])/2),-(((Sqrt[1-(cos^2)]) (Sqrt[s (-4 (M^2)+s)]) Cos[ϕ])/(2 (Sqrt[s]))),-(((Sqrt[1-(cos^2)]) (Sqrt[s (-4 (M^2)+s)]) Sin[ϕ])/(2 (Sqrt[s]))),-((cos (Sqrt[s (-4 (M^2)+s)]))/(2 (Sqrt[s])))}
};

We can verify that s=((((Subscript[p, 1])+(Subscript[p, 2])))^2), (m^2)=(Subscript[p, 1]^2)=(Subscript[p, 2]^2), and (M^2)=(Subscript[p, 3]^2)=(Subscript[p, 4]^2)

{
s==sq[p1+p2]==sq[p3+p4],
(m^2)==sq[p1]==sq[p2],
(M^2)==sq[p3]==sq[p4]
}/.sq[p_]:>p.p/.momenta//Simplify

More importantly, we can calculate t=((((Subscript[p, 1])-(Subscript[p, 3])))^2)

tsol=Simplify[t->sq[p1-p3]/.sq[p_]:>p.p/.momenta]

Our matrix element can be written as

ℳsqCos=Simplify[ℳsqS/.momenta]

The velocity of the muon is defined as

velcond=β==Hold[((Norm[p3[[2;;4]]])/(p3[[1]]))]/.momenta//ReleaseHold//Simplify

Step 11-12: cross section

We can now calculate the cross section as

calculate the cross section as

dσ=(1/((2(Subscript[E, 1]))(2(Subscript[E, 2]))(Abs[(OverVector[v])])))(Underscript[(Underscript[((∏_(n=3)^f)(((d^3)(Subscript[p, n]))/((((2π))^3)2(Subscript[E, n]))))(((2π))^4)(δ^((4)))((Subscript[p, 1])+(Subscript[p, 2])-(∑_(n=3)^f)(Subscript[p, i])), ⏟]), d(Subscript[Φ, 2])])(1/4)∑((Abs[ℳ])^2)

Using what we know about the phase space

d(Subscript[Φ, 2])=dΩ (1/(16(π^2)))((Abs[(OverVector[(Subscript[p, 3])])])/(Subscript[E, cm]))

We have

(1/((2E1)(2E2)v))dΦ (1/4)ℳsqCos/.dΦ->(1/(16(π^2)))((Norm[p3])/Ecm)
dσdΩ=%//.{
(1/(v E1 E2))->(1/(Sqrt[(((p1.p2))^2)-(m^2)(m^2)])),
Ecm->E1+E2,
E1->Hold[p1[[1]]],
E2->Hold[p2[[1]]],
Norm[p3]->Norm[Hold[p3[[2;;4]]]]
}/.momenta
dσdΩ=Simplify[ReleaseHold[%]/.e->(Sqrt[4π α])]

In the lecture we then set m=0 and identified the velocity of the muon β

dσdΩm0=Simplify[dσdΩ/.Last@Solve[velcond,M]/.m->0]

We can now set for example the muon mass to zero, i.e. β=1

dσdΩm0M0=Simplify[dσdΩ/.m|M->0]

or calculate the total cross section

integrate over phase space

σtot=(∫_0^(2π))(∫_(-1)^1)dσdΩ\[DifferentialD]cos\[DifferentialD]ϕ

Plots

We can now show our result. For example at the KLOE experiment which has (Sqrt[s])=1.02 GeV, we have

Quantity[1,"SpeedOfLight""ReducedPlanckConstant"]

numbers={
conv->((Quantity[1,"SpeedOfLight""ReducedPlanckConstant"])^2),
s->((Quantity[1.02,"Gigaelectronvolts"])^2),
m->Quantity[0.51099895,"Megaelectronvolts"],
M->Quantity[105.658375,"Megaelectronvolts"],
α->(1/137.035999084),
cos->Quantity[cos,1]
};

dσdΩSI=Simplify@UnitConvert[Simplify[conv dσdΩ/.numbers],"Nanobarns"]

Plot[dσdΩSI,{cos,-1,1},PlotRange->All]

UnitConvert[conv σtot/.numbers,"Nanobarns"]

Going beyond the lecture (not examinable)

We can also calculate the one-loop corrections to this.

amplitude1=RunQGraf[{"elm","elp"},{"mum","mup"},1][[2;;,2]]/.line[a_,b_]:>b/._δZ|_δm->0

GraphicsGrid@Partition[QGrafGetDiagrams["qgraf[{elm,elp},{mum,mup},1]"][[2;;]],2]

Note that we now have both closed fermion loops (first row of diagrams) and unconstrained momenta. The generator chose a specific routing for us, i.e. it identified a momentum to be the loop momentum (Subscript[k, 1]). We can use the technology we are developing in Chapter 7 to calculate this in dimensional regularisation

Export for the web

Export[NotebookDirectory[]<>"/ee2mm.html", InputNotebook[]]