QED calculation: 1-Loop Calculation

Initialization

Quit[];

$Language="English";
<<FeynArts`
<<FormCalc`

(* helicity to be averaged *)
_Hel=0;

Create Topology

topo0=CreateTopologies[0,2->2];
topo1=CreateTopologies[1,2->2,ExcludeTopologies->{Internal}];
topo01=Join[topo0,topo1];
Paint[topo01, ColumnsXRows->4];

(SuperPlus[e])(SuperMinus[e])->(SuperPlus[μ])(SuperMinus[μ]): QED tree-level and1-loop diagrams

process={F[2,{1}],-F[2,{1}]}->{F[2,{2}],-F[2,{2}]};
TreeDiagrams = InsertFields[topo0,process,InsertionLevel->{Classes},
ExcludeParticles->{V[2],V[3],V[4],S}];
OneLoopDiagrams = InsertFields[topo1,process,InsertionLevel->{Classes},
ExcludeParticles->{V[2],V[3],V[4],S}];
FullDiagrams=Join[TreeDiagrams,OneLoopDiagrams];
Paint[FullDiagrams];

Tree-level amplitude

ClearProcess[];
TreeMatrixSquared=Module[{amp,sq,mat},
amp = CalcFeynAmp[CreateFeynAmp[TreeDiagrams],FermionChains->VA];
sq =SquaredME[amp]//. HelicityME[amp];
mat =  sq[[1]]//.sq[[2]];
FullSimplify[4*mat /. Den[a_, b_] :> (1/(a-b)),S+T+U==2ME2+2MM2]]

Then we will calculate one-loop level amplitude, but notice that we have to execute “ClearProcess[]” before the one-loop level evaluation.
It will clear all the abbreviations / subexpressions, so we have to be sure that our tree-level result does not contain abbreviations and subexpressions.

Here, our result, which is stored as TreeMatrixSquared, has no abbreviations, so we are ready to execute ClearProcess[].

Tree + 1-loop level amplitude

ClearProcess[];
FullAmp = CalcFeynAmp[CreateFeynAmp[FullDiagrams],FermionChains->VA]

FullSq = SquaredME[FullAmp]//.HelicityME[FullAmp];
FullMatrixSquared =4*FullSq[[1]]//.FullSq[[2]]

Numerical Evaluation with LoopTools

You can load LoopTools by:

Install["LoopTools"]

Now we can evaluate the Passarino-Veltman integrals with LoopTools.
For example,

D0i[dd1,5,1,-2,1,1,1,1,1,1,1]

So,

replaces=Join[
Subexpr[],
{Den[a_,b_]:>1/(a-b)},
{S->(((2EE))^2),T->-(2(EE^2)-ME2-MM2)+2(Sqrt[(EE^2)-ME2])(Sqrt[(EE^2)-MM2])Cosθ,U->2ME2+2MM2-S-T},
{EE->125,ME->0.511*^-3,ME2->(((0.511*^-3))^2),MM->0.106,MM2->(0.106^2),Alfa->α,Alfa2->(α^2)}];
{TreeMatrixSquared,FullMatrixSquared}//.replaces/.Cosθ->0.5//FullSimplify

Note that the IR-divergence is detected.
You should be aware of the IR/UV divergences BEFORE you write computation code!

We can verify the presence of IR divergence by changing (λ^2) (fictitious photon mass).
The result does depend on (λ^2), so IR divergence is not cancelled.

SetLambda[1];
Chop[FullSimplify[FullMatrixSquared//.replaces/.Cosθ->0.5]]
SetLambda[(0.1^2)];
Chop[FullSimplify[FullMatrixSquared//.replaces/.Cosθ->0.5]]
SetLambda[(10^2)];
Chop[FullSimplify[FullMatrixSquared//.replaces/.Cosθ->0.5]]

In dimensional regularization, it will be

constant=FullMatrixSquared/.{C0i[__]->0,D0i[__]->0}//.replaces/.Cosθ->0.5;

SetLambda[-2];
pole2=FullMatrixSquared//.replaces/.Cosθ->0.5;
SetLambda[-1];
pole1=FullMatrixSquared//.replaces/.Cosθ->0.5;
SetLambda[0];
finite=FullMatrixSquared//.replaces/.Cosθ->0.5;
Collect[Chop[Expand[((pole2-constant)/(ϵ^2))+((pole1-constant)/ϵ)+finite]],α]//Expand