Sample FeynCalc notebook

Examples of  Mathematica calculations of  amplitudes, cross sections, and identities in the text MouseOver using the FeynCalc package. (Tested using FeynCalc vs 9.0.1 with Mathematica 10.4 and 11.0 on Macintosh.) Caveat: there is no guarantee that the commands used are the most elegant or most efficient. (Simple calculations not requiring FeynCalc are illustrated in matrixelements.nb.)

FeynCalc writeups at http://arxiv.org/pdf/1604.06709.pdf and http://arxiv.org/pdf/1601.01167.pdf. Download from
https://github.com/FeynCalc/feyncalc/releases or uncomment and run the following:

(*Import["https://raw.githubusercontent.com/FeynCalc/feyncalc/master/install.m"]
InstallFeynCalc[]*)

Metric and gamma matrix conventions are as in Bjorken and Drell. However, (ϵ^0123)=+1 is opposite of BD (and my) convention. Spinors are covariantly normalized. Gauge couplings have opposite of SMB sign. Most calculations in this notebook use internally defined vertices, etc, rather than built-in ones, to conform with SMB conventions.

May need to  reload FeynCalc (or even to first delete all output) to clear variables between calculations. Also, sometimes helps to restart kernel.

The Initialize section loads FeynCalc, defines macros, sets options, etc. Ignore error messages.

Go to Macros, Traces, QED, QCD, Fermi Theory, EW Model

Clear["Global`*"]
DateString[]
{$Version,$UserName}
<<FeynCalc`

Macros

Define general field theory and QCD parameters and macros

SUNN=3; (*set number of colors*)
SetOptions[SUNSimplify,SUNNToCACF->False];
(*{SUNSimplify[CF],SUNSimplify[CA]}*)
GlPr[p_,μ_,i_,ν_,j_]:=-I  MT[μ,ν] SUNDelta[i,j]/ScalarProduct[p]//ExpandScalarProduct
(*GlPr[p_,μ_,i_,ν_,j_]:=GluonPropagator[Momentum[p],μ,i,ν,j,Gauge->1,Dimension->4]//ExpandScalarProduct//Explicit*)
ScPr[p_]:=I/ScalarProduct[p]//ExpandScalarProduct (*Scalar Propagator. Use -ScPr[p] for vector without indices*)
ScPrm[p_,m_]:=I/(ScalarProduct[p]-m^2)//ExpandScalarProduct 
QPr[p_,m_]:= I (DiracSlash[p]+m)/(ScalarProduct[p]-(m^2))//ExpandScalarProduct
(* Introduce since GluonPropagator and QuarkPropagator do not substitute in denominators. One could also use PropagatorDenominatorExplicit for this purpose*)
cgv[p_,μ_,q_,ν_,r_,σ_]:=MT[μ,ν]FourVector[q-p,σ]+MT[μ,σ]FourVector[p-r,ν]+MT[ν,σ]FourVector[r-q,μ] (*triple gauge vertex without group indices. Full QCD vertex (with Feyncalc sign convention) is
 -(Subscript[c, ijk]) Gstrong cgv[],
which is equivalent to GluonVertex[]*)
qgv[μ_,i_]:=-QGV[μ,i] (*convert to SMB sign convention (optional) *)

{qav=6,Gav=16, eav=2};
PrSpin[stn_]:=(1+ GA[5].GS[stn])/2   (*spin projection. Use GA[7] for massless (Subscript[P, L]), GA[6] for massless (Subscript[P, R])*)

(*----------version 9 additions and changes---------*)
(*ftrace ={DiracTrace->Tr2,D->4};*) (*use with color calculations*)
ftrace =DiracTrace->TR;  (*modify for vs 9*)
(*Replace QGV[μ,i] by QGV9[μ,i,a,b] throughout*)
(Subscript[g, s])= Gstrong;   (*this ensures that explicit (Subscript[g, s])is recongnized as Gstrong*)
QGV9[μ_,i_,a_,b_]:= -I (Subscript[g, s])GA[μ]SUNTF[{i},a,b]  (*new for vs 9; handlesquark color indices correctly and uses SMB sign convention *)
QGV9[μ,i,a,b];
(*Change ClearScalarProducts to FCClearScalarProducts[]*)
GV9[{p_,μ_,i_},{q_,ν_,j_},{r_,σ_,k_}]:=(Subscript[g, s])  SUNF[i,j,k]( MT[μ,ν] FourVector[q-p,σ]+ MT[μ,σ] FourVector[p-r,ν]+ MT[ν,σ] FourVector[r-q,μ])   (*triple gluon vertex using SMB conventions*)
GV9[{p,μ,i},{q,ν,j},{r,σ,k}];
GV49[{μ_,i_},{ν_,j_},{ρ_,l_},{σ_,k_}]:= -I (Subscript[g, s])^2 (SUNF[i,j,m]SUNF[k,l,m](MT[μ,σ]MT[ν,ρ]-MT[μ,ρ]MT[ν,σ] )+
SUNF[i,k,m]SUNF[j,l,m](MT[μ,ν]MT[ρ,σ]-MT[μ,ρ]MT[ν,σ] )+
SUNF[i,l,m]SUNF[j,k,m](MT[μ,ν]MT[ρ,σ]-MT[μ,σ]MT[ν,ρ] ))
(*quartic gauge vertex*)
GV49[{μ,i},{ν,j},{ρ,l},{σ,k}];
GV49a[{p1_,μ_,i_},{p2_,ν_,j_},{p3_,ρ_,l_},{p4_,σ_,k_}]:= GV49[{μ,i},{ν,j},{ρ,l},{σ,k}]   
GV49a[{p1,μ,i},{p2,ν,j},{p3,ρ,l},{p4,σ,k}];(*include irrelevant momenta*)
(*QsGV[p_,μ_,i_]:=-Gstrong I FourVector[p,μ]SUNT[i]*)  (*scalar gluon vertex in SMB sign convention*)(*momentum is (Subscript[(OverBar[p]), a])+(Subscript[(OverBar[p]), b]) figure 4.1, p162*)

QsGV9[p_,μ_,i_,a_,b_]:=-Gstrong I FourVector[p,μ]SUNTF[{i},a,b];  (*new for vs 9*)
QsGV9[p,μ,i,a,b];
(*QsQsGGV[{μ_,i_},{ν_,j_}]:=I Gstrong^2 (SUNT[i].SUNT[j]+SUNT[j].SUNT[i]) MT[μ,ν]*)(*GG (Subscript[q, 0]) (Subscript[q, 0]) vertex*)

QsQsGGV9[{μ_,i_},{ν_,j_},a_,b_]:=I Gstrong^2 (SUNTF[{i,j},a,b]+SUNTF[{j,i},a,b]) MT[μ,ν];
QsQsGGV9[{μ,i},{ν,j},a,b];

Define SM parameters and macros

Use parameters e,  g , g', (Subscript[g, Z]),(Subscript[s, W]), (Subscript[c, W]), (Subscript[M, W]), (Subscript[M, Z]), (Subscript[M, H]), (Subscript[G, F]), etc, as in SMB.
Will work with SMB sign convention for gauge couplings.

Define SM vertices

EwPr[p_,M_]:=-I/(ScalarProduct[p]-M^2)//ExpandScalarProduct 
(*aww is equivalent to SMVertex["AWW", p, μ, q, ν, r, σ] except for the sign convention.
Could also use built in SMVertex["HHH",__] *)
(*test*)(*{aww[p, μ, q, ν, r, σ],zww[p, μ, q, ν, r, σ],wfbfL[μ],wfbfR[μ],wfbf[μ],afbfL[μ,-1],afbfR[μ,-1],afbf[μ,-1],zfbfL[μ,-1/2,-1],zfbfR[μ,-1/2,-1],zfbf[μ,-1/2,-1],EwPr[p1+p2,(Subscript[M, Z])]}*)
BWPr[p_,M_,Γ_]:=-I/(ScalarProduct[p]-M^2+ I M Γ)//ExpandScalarProduct 
aww[p_, μ_, q_, ν_, r_, σ_]:=+ I g sw cgv[p, μ, q, ν, r, σ]
zww[p_, μ_, q_, ν_, r_, σ_]:=+ I g cw cgv[p, μ, q, ν, r, σ]
zzH[μ_, ν_]:=+ I g MZ MT[μ,ν]/cw 
wwH[μ_, ν_]:=+ I g MW MT[μ,ν]
fbfH[mfv_]:=-I mfv/nuH
HHH[mhv_]:=-3 I mhv^2/nuH
wfbfL[μ_]:=- I g GA[μ].GA[7] /Sqrt[2]
wfbfR[μ_]:=0
wfbf[μ_]:=wfbfL[μ]+wfbfR[μ]
afbfL[μ_,qv_]:=- I g  sw qv GA[μ].GA[7]
afbfR[μ_,qv_]:=- I g sw qv GA[μ].GA[6]
afbf[μ_,qv_]:=- I g sw qv GA[μ]
γfbf[μ_,qv_]:=- I Sqrt[4 Pi α]  qv GA[μ]    (*use for pure QED *)
(* Z parameters and couplings*)
parv={t3lv={1/2,-1/2,1/2,-1/2},qvv={0,-1,2/3,-1/3},COLv={1,1,3,3}};
gVv[i_]:=t3lv[[i]]-2 sw^2 qvv[[i]];gAv[i_]:=t3lv[[i]]; 
epslv[i_]:=t3lv[[i]]- sw^2 qvv[[i]]; epsrv[i_]:=- sw^2 qvv[[i]](*1=ν,2=e,3=u,4=d*)
zfbfL[μ_,t3l_,qv_]:=- I g  (t3l-sw^2 qv)GA[μ].GA[7]/cw  
zfbfR[μ_,t3l_,qv_]:=- I g  (-sw^2 qv)GA[μ].GA[6]/cw
zfbf[μ_,t3l_,qv_]:=zfbfL[μ,t3l,qv]+zfbfR[μ,t3l,qv]
zfbfg[μ_,gV_,gA_]:=-I g  GA[μ].(gV-gA  GA[5])/(2 cw)
zfbfϵ[μ_,epsl_,epsr_]:=-I g  GA[μ].(epsl GA[7]+ epsr  GA[6])/ cw
intofin={gV->gVf,gA->gAf,epsl->epslf,epsr->epsrf,t3l->t3lf, qv->qvf, COL->COLf};
intok[k_]:={gV->gVv[k],gA->gAv[k],epsl->epslv[k],epsr->epsrv[k],t3l->t3lv[[k]], qv->qvv[[k]], COL->COLv[[k]]}
fintok[k_]:={gVf->gVv[k],gAf->gAv[k],epslf->epslv[k],epsrf->epsrv[k],t3lf->t3lv[[k]], qvf->qvv[[k]], COLf->COLv[[k]]}
ncpar={epsl-> t3l-sw^2 qv, epsr->-sw^2 qv, gV-> t3l-2 sw^2 qv,gA-> t3l,epslf-> t3lf-sw^2 qvf, epsrf->-sw^2 qvf, gVf-> t3lf-2 sw^2 qvf,gAf-> t3lf};
asympar={gAf gVf-> Af (gAf^2+gVf^2)/2 ,gA gV-> A (gA^2+gV^2)/2};
display = {MZ->(Subscript[M, Z]),GZ->(Subscript[Γ, Z]), MW->(Subscript[M, W]), mt->(Subscript[m, t]),GF-> (Subscript[G, F]),MH->(Subscript[M, H]),nuH->ν,sw->(Subscript[s, W]),cw->(Subscript[c, W]),tw->(Subscript[t, W]),sz->(Subscript[(OverHat[s]), Z]),cz->(Subscript[(OverHat[c]), Z]),ghat->(OverHat[g]),gV->(Subscript[g, V]),gA->(Subscript[g, A]),epsl->(Subscript[ϵ, L]),epsr->(Subscript[ϵ, R]),t3l->(Subscript[t, L]^3), qv->q, COL->C3,gVf->(Subscript[g, Vf]),gAf->(Subscript[g, Af]),epslf->(Subscript[ϵ, L]^f),epsrf->(Subscript[ϵ, R]^f),t3lf->(Subscript[t, fL]^3), qvf->(Subscript[q, f]), COLf->(Subscript[C3, f]),Af->(Subscript[A, f]),cθ->(Subscript[c, θ])};
weakangle={sw->Sin[(Subscript[θ, W])], cw->Cos[(Subscript[θ, W])],tw->Tan[(Subscript[θ, W])]};
smpar={ (Subscript[G, F])->1.17 10^(-5), GF->1.17 10^(-5),α->1/137,(Subscript[α, W])->1/128,g->.654, (Subscript[θ, W])-> .501,cw-> 0.88,sw->0.48,tw->0.55, MZ->91.19, GZ->2.50,MW->80.4, (Subscript[M, W])->80.4, (Subscript[M, Z])-> 91.19, hbcsq->(.197)^2 10^(10) (*pb-GeV^2*),hbar->6.6 10^(-25) (*GeV-s*),hbarc->0.197 10^(-13) (*GeV-cm*),clight->3.0 10^(10) (*cm/s*), mt->173.1,nuH->246};

Kinematics macros

kfc[s_,M1_,M2_]:=Sqrt[s-(M1+M2)^2]Sqrt[s-(M1-M2)^2]/(2 Sqrt[s])
e1fc[s_,M1_,M2_]:=(s+M1^2-M2^2)/(2 Sqrt[s])

2 x 2 scattering

Following command defines general kinematics for 2 x 2 scattering.

(*FCClearScalarProducts[];
SetMandelstam[s, t, u, pg1,pg2, -pg3, -pg4,mg1,mg2,mg3,mg4];
{kgf=kfc[s,mg3,mg4],kgi=kfc[s,mg1,mg2],eg1=e1fc[s,mg1,mg2],eg2=e1fc[s,mg2,mg1],eg3=e1fc[s,mg3,mg4],eg4=e1fc[s,mg4,mg3],βg1=kgi/eg1,βg2=kgi/eg2, βg3=kgf/eg3, βg4=kgf/eg4};
coefg=kgf/(32 Pi s kgi); (*coefficient for ((dσ)/dz) *)
cmvar={t->mg1^2+mg3^2-2 eg1 eg3 + 2 kgi kgf cθ, u->mg1^2+mg4^2-2 eg1 eg4 - 2 kgi kgf cθ};
{thresholdi={s->(mg1+mg2)^2},thresholdf={s->(mg3+mg4)^2}};
ctheta = {cθ->Cos[(Subscript[θ, CM])]};
*)

Can follow by redefining masses, e.g.,

(*{mg1=M,mg2=M,mg3=MZ,mg4=MH}*)

2 body decays

Following command defines general kinematics for 2 body decays. Copy to start of calculation.

(*
FCClearScalarProducts[];
{ScalarProduct[pdl]=md1^2,ScalarProduct[pd2]=md2^2,ScalarProduct[pd3]=md3^2,ScalarProduct[pd2,pd3]=(md1^2-md2^2-md3^2)/2,ScalarProduct[pd1,pd3]=(md1^2+3md3^2-md2^2)/2,ScalarProduct[pd1,pd2]=(md1^2+3md2^2-md2^3)/2};
{kdf=kfc[md1^2,md2,md3],ed2=e1fc[md1^2,md2,md3],ed3=e1fc[md1^2,md3,md2],βd2=kdf/ed2,βd3=kdf/ed3}//Refine[#,md1>0]&;
coefd=kdf/(8 Pi md1^2);  (*coefficient for total Γ *)
{ScalarProduct[sd1]=-1,ScalarProduct[sd1,pd1]=0,ScalarProduct[sd1,pd2]=-kdf cθ,ScalarProduct[sd1,pd3]=+kdf cθ};
ctheta = {cθ->Cos[(Subscript[θ, CM])]};
*)

Above sets up general kinematics for 2 body decays. Can follow by redefining masses, e.g.,

(*{md1=M,md2=MZ,md3=MZ}*)

TOC

Trace Identities, p42

TR[GS[a].GS[b]]
TR[GS[a].GS[b].GS[c].GS[d]]
TR[GS[a].GS[b].GS[c].GS[d].GS[f].GS[g]] 
TR[GS[a].GS[b].GS[c].GS[d].GA[5]]
SetOptions[DiracTrace,LeviCivitaSign->+1 ];(*default is -1*)
TR[GS[a].GS[b].GS[c].GS[d].GA[5]]
SetOptions[DiracTrace,LeviCivitaSign->-1 ];
GA[μ,ν,ρ]
Chisholm[%]
EpsChisholm[%]
TR[GS[a].GS[b].GS[c].GS[d].GS[f].GS[g].GA[5]]

TR[GA[μ].GS[a].GA[ν].GS[b].(1+GA[5])]TR[GA[μ].GS[c].GA[ν].GS[d].(1+λ GA[5])]//Contract//Simplify
%/. λ->+1 //Simplify
%%/. λ->-1//Simplify

{GA[μ].GA[μ]//DiracSimplify,GA[μ].GS[a].GA[μ]//DiracSimplify,GA[μ].GS[a].GS[b].GA[μ]//DiracReduce,GA[μ].GS[a].GS[b].GA[c].GA[μ]//DiracSimplify,GA[μ].GS[a].GS[b].GA[c].GA[d].GA[μ]//DiracSimplify}

LeviCivita[μ,ν,ρ,σ]LeviCivita[μ,ν,ρ,σ]//Contract
LeviCivita[μ,ν,ρ,σ]LeviCivita[μ1,ν,ρ,σ]//Contract
LeviCivita[μ,ν,ρ,σ]LeviCivita[μ1,ν1,ρ,σ]//Contract
LeviCivita[μ,ν,ρ,σ]LeviCivita[μ1,ν1,ρ1,σ]//Contract
LeviCivita[μ,ν,ρ][a]LeviCivita[μ,ν,ρ][b]//Contract
LeviCivita[μ,ν][a,b]LeviCivita[μ,ν][c,d]//Contract

TOC

QED Calculations

(SuperMinus[e])(SuperPlus[e]) → f (OverBar[f]) from p 53

Define kinematics

FCClearScalarProducts[];
{spcase={(Subscript[m, i])->0, (Subscript[m, f]) -> Sqrt[1-((Subscript[β, f])^2)](Subscript[E, cm])}, cme={s-> 4Ε^2}, plotpar={α->1/137,(Subscript[Q, f])->1}};  (*special cases*)
{(Subscript[E, cm])=Sqrt[s]/2, (Subscript[p, i])=Sqrt[(Subscript[E, cm])^2-(Subscript[m, i])^2], (Subscript[p, f])=Sqrt[(Subscript[E, cm])^2-(Subscript[m, f])^2],
coef0=-I 4 Pi α  (Subscript[Q, f])/s, coef1=4 ((Subscript[Q, f])^2)(Pi^2) (α^2)/(s^2), coef2= (Subscript[p, f])/(32  Pi s (Subscript[p, i]))};
{coef3=coef2 /. spcase// Cancel ,coef4=Refine[coef3,{s>0, (Subscript[β, f])>0}]};
{ScalarProduct[p1,p1]=ScalarProduct[p2,p2]=((Subscript[m, i])^2),
ScalarProduct[p3,p3]=ScalarProduct[p4,p4]=((Subscript[m, f])^2),
ScalarProduct[p1,p2]=2(Subscript[E, cm])^2-(Subscript[m, i])^2,ScalarProduct[p3,p4]=2(Subscript[E, cm])^2-(Subscript[m, f])^2,
ScalarProduct[p1,p4]=ScalarProduct[p2,p3]= ((Subscript[E, cm])^2+(Subscript[p, i])  (Subscript[p, f])  z),
ScalarProduct[p1,p3]=ScalarProduct[p2,p4]= ((Subscript[E, cm])^2-(Subscript[p, i])  (Subscript[p, f])  z)};
{{p1sl,p2sl,p3sl,p4sl}=DiracSlash/@{p1,p2,p3,p4}};

Square the amplitude by hand (Included for illustration. Not needed in calculation.)

tra=TR[GA[ν].(p4sl-(Subscript[m, f])).GA[μ].(p3sl+(Subscript[m, f]))] TR[GA[ν].(p1sl+(Subscript[m, i])).GA[μ].(p2sl-(Subscript[m, i]))]
prod=Simplify[Contract[tra]/.spcase]
coef1 prod

Now work with amplitude directly instead

(Subscript[M, fi])=coef0 SpinorUBar[p3,(Subscript[m, f])].GA[μ].SpinorV[p4,(Subscript[m, f])] SpinorVBar[p2,(Subscript[m, i])].GA[μ].SpinorU[p1,(Subscript[m, i])]//FCI
(Subscript[Mstar, fi])=ComplexConjugate[(Subscript[M, fi])]/.μ->ν
(Subscript[Mbarsq, fi])=(1/4) FermionSpinSum[(Subscript[M, fi]) (Subscript[Mstar, fi])]
Contract[%/.DiracTrace->TR]
Mbarsq=Simplify[%/.spcase]

Calculate cross section

dsigbardz=coef4 Mbarsq
sigbar=Integrate[dsigbardz,{z,-1,1}]//Simplify
{dsigbarreldz=dsigbardz,sigbarrel=sigbar}/.(Subscript[β, f])->1

shape=((100 Mbarsq)/(Subscript[β, f]))/.plotpar;
Plot[{shape/.(Subscript[β, f])->0.5,shape/.(Subscript[β, f])->0.75,shape/.(Subscript[β, f])->1},{z,-1,1}]

Bhabha Scattering, p54

ignore masses

FCClearScalarProducts[];
{coef0=4 Pi I α, coef1=1/(32 Pi s)};
{ScalarProduct[p1a,p1a]=ScalarProduct[p2a,p2a]=0,
ScalarProduct[p3a,p3]=ScalarProduct[p4a,p4a]=0,
ScalarProduct[p1a,p2a]=ScalarProduct[p3a,p4a]=s/2,ScalarProduct[p1a,p4a]=ScalarProduct[p2a,p3a]= -u/2,
ScalarProduct[p1a,p3a]=ScalarProduct[p2a,p4a]= -t/2};
(* {{p1sl,p2sl,p3sl,p4sl}=DiracSlash/@{p1,p2,p3,p4}}; *)

f1=coef0 SpinorUBar[p3a,0].GA[μ].SpinorV[p4a,0] SpinorVBar[p2a,0].GA[μ]. SpinorU[p1a,0] //FCI
f2=coef0 SpinorUBar[p3a,0].GA[μ].SpinorU[p1a,0] SpinorVBar[p2a,0].GA[μ]. SpinorV[p4a,0] //FCI
(* Could also combine the two amplitudes at this stage, as in next example *)
f1star=ComplexConjugate[f1] /. μ-> ν 
f2star=ComplexConjugate[f2] /. μ-> ν 
m11= FermionSpinSum[f1 f1star] /.DiracTrace->TR// Contract //Simplify;
m22= FermionSpinSum[f2 f2star] /.DiracTrace->TR// Contract //Simplify;
m12= FermionSpinSum[f1 f2star] /.DiracTrace->TR// Contract //Simplify;
Mbarsq= (m11/(s^2)+m22/(t^2)-2  m12/(s t))/4  //Factor (*uses reality*);
dsigbardz=coef1  Mbarsq ;
dsigbardz/(Pi (α^2)/s) //Expand;
Collect[C2 %,C2] /. C2->(Pi (α^2)/s)
(*TrickMandelstam[%,{s,t,u,0}]*)

Repeat Bhabha scattering, keeping common mass m and using SetMandelstam

Keep common mass m. Notes: (a) SetMandelstam assumes all momenta are incoming; (b) to clear the SetMandelstam assignments must delete all ouput and restart. We avoid that here by using different names. Carry out calculation in two forms.

FCClearScalarProducts[];
{coef0=4 Pi I α, coef1=1/(32 Pi s)};SetMandelstam[s, t, u, p1b, p2b, -p3b, -p4b,m,m,m,m];
(* {{p1bsl,p2bsl,p3bsl,p4bsl}=DiracSlash/@{p1b,p2b,p3b,p4b}}; *)

Next command handles the two terms in amplitude separately

f1=coef0 SpinorUBar[p3b,m].GA[μ].SpinorV[p4b,m] SpinorVBar[p2b,m].GA[μ]. SpinorU[p1b,m] //FCI
f2=coef0 SpinorUBar[p3b,m].GA[μ].SpinorU[p1b,m] SpinorVBar[p2b,m].GA[μ]. SpinorV[p4b,m] //FCI
f1star=ComplexConjugate[f1] /. μ-> ν ;
f2star=ComplexConjugate[f2] /. μ-> ν ;
m11= FermionSpinSum[f1 f1star] /.DiracTrace->TR// Contract //Simplify;
m22= FermionSpinSum[f2 f2star] /.DiracTrace->TR// Contract //Simplify;
m12= FermionSpinSum[f1 f2star] /.DiracTrace->TR// Contract //Simplify;
Mbarsq= (m11/(s^2)+m22/(t^2)-2  m12/(s t))/4  //Factor; (*uses reality*)
dsigbardz=coef1  Mbarsq ;
dsigbardz/(Pi (α^2)/s) //Expand;
Collect[C2 %,C2] /. C2->(Pi (α^2)/s)
 TrickMandelstam[%,{s,t,u,4 (m^2)}]
%/.m->0

Simpler form combines the two amplitudes

amp= coef0 ( SpinorUBar[p3b,m].GA[μ].SpinorV[p4b,m] SpinorVBar[p2b,m].GA[μ]. SpinorU[p1b,m]/s-  SpinorUBar[p3b,m].GA[μ].SpinorU[p1b,m] SpinorVBar[p2b,m].GA[μ]. SpinorV[p4b,m]/t)//FCI
ampst= ComplexConjugate[amp] /. μ-> ν ;
Mbsq =FermionSpinSum[amp ampst/4 //Expand] /.DiracTrace->TR // Contract //Simplify //Factor;
  (*need to Expand amp ampst before FermionSpinSum*)
dsigbdz=coef1  Mbsq ;
dsigbdz/(Pi (α^2)/s) //Expand;
Collect[C3 %,C3] /. C3->(Pi (α^2)/s)
%/. m->0

Moller scattering, keeping common mass m and using SetMandelstam

Keep common mass m. Notes: (a) SetMandelstam assumes all momenta are incoming; (b) to clear the SetMandelstam assignments must delete all ouput and restart. We avoid that here by using different names. Carry out calculation in two forms.

FCClearScalarProducts[];
{coef0=4 Pi I α, coef1=1/(32 Pi s)};SetMandelstam[s, t, u, p1b, p2b, -p3b, -p4b,m,m,m,m];
(* {{p1bsl,p2bsl,p3bsl,p4bsl}=DiracSlash/@{p1b,p2b,p3b,p4b}}; *)

Next command handles the two terms in amplitude separately

f1=coef0 SpinorUBar[p3b,m].GA[μ].SpinorU[p1b,m] SpinorUBar[p4b,m].GA[μ]. SpinorU[p2b,m] //FCI
f2=coef0 SpinorUBar[p4b,m].GA[μ].SpinorU[p1b,m] SpinorUBar[p3b,m].GA[μ]. SpinorU[p2b,m] //FCI
f1star=ComplexConjugate[f1] /. μ-> ν ;
f2star=ComplexConjugate[f2] /. μ-> ν ;
m11= FermionSpinSum[f1 f1star] /.DiracTrace->TR// Contract //Simplify;
m22= FermionSpinSum[f2 f2star] /.DiracTrace->TR// Contract //Simplify;
m12= FermionSpinSum[f1 f2star] /.DiracTrace->TR// Contract //Simplify;
Mbarsq= (m11/(t^2)+m22/(u^2)-2  m12/(u t))/4  //Factor; (*uses reality*)
dsigbardz=coef1  Mbarsq ;
dsigbardz/(Pi (α^2)/s) //Expand;
Collect[C2 %,C2] /. C2->(Pi (α^2)/s)
 TrickMandelstam[%,{s,t,u,4 (m^2)}]
%/.m->0

Simpler form combines the two amplitudes

amp= coef0 (SpinorUBar[p3b,m].GA[μ].SpinorU[p1b,m] SpinorUBar[p4b,m].GA[μ]. SpinorU[p2b,m] /t-   SpinorUBar[p4b,m].GA[μ].SpinorU[p1b,m] SpinorUBar[p3b,m].GA[μ]. SpinorU[p2b,m]/u)//FCI
ampst= ComplexConjugate[amp] /. μ-> ν ;
Mbsq =FermionSpinSum[amp ampst/4 //Expand] /.DiracTrace->TR // Contract //Simplify //Factor;
  (*need to Expand amp ampst before FermionSpinSum*)
dsigbdz=coef1  Mbsq ;
dsigbdz/(Pi (α^2)/s) //Expand;
Collect[C3 %,C3] /. C3->(Pi (α^2)/s)
%/. m->0

Pion Compton Scattering, p38

Define polarizations in lab frame

FCClearScalarProducts[];
coef0=-I 4 Pi α; coef1=(k2^2)/(32 Pi (m^2)(k1^2));
{ep1p=Polarization[k1f,Transversality->True],
ep1=Momentum[ep1p],
ep2p=Polarization[k2f,Transversality->True],
ep2=Momentum[ep2p],
ep1ps= ComplexConjugate[ep1p],
ep1s=Momentum[ep1ps],
ep2ps= ComplexConjugate[ep2p],
ep2s=Momentum[ep2ps]};
{{p1v,p2v,k1v,k2v}=Momentum/@{p1f,p2f,k1f,k2f}};
{ScalarProduct[p1f,p1f]=ScalarProduct[p2f,p2f]=(m^2),
ScalarProduct[k1f,k1f]=ScalarProduct[k2f,k2f]=0,
ScalarProduct[k1f,k2f]=k1 k2 (1-zl),ScalarProduct[p1f,k1vf]=m k1,ScalarProduct[p1f,k2f]=m k2,
ScalarProduct[p1f,p2f]=m (m+k1-k2), ScalarProduct[k1f,p2f]=m k1 - k1 k2 (1-zl), ScalarProduct[k2f,p2f]=m k2 + k1 k2 (1-zl),
ScalarProduct[p1f,ep1p]=ScalarProduct[p1f,ep2p]=0,
ScalarProduct[p1f,ep1ps]=ScalarProduct[p1f,ep2ps]=0};
pair1=(Pair[ep2s,k2v]+2 Pair[ep2s,p2v])(Pair[ep1,k1v]+2Pair[ep1,p1v])
pair2=(2Pair[ep2s,p1v]-Pair[ep2s,k2v])(2 Pair[ep1,p2v]-Pair[ep1,k1v]);
pair3=-2 Pair[ep2s,ep1];
coef0 (pair1/(2 m k1)- pair2/(2 m k2) + pair3)
(Subscript[M, fi])= %/. Pair[ep1,ep2s]->S

Differential cross section in lab for massless photon

coef1 (Subscript[M, fi]) (Subscript[M, fi])\[Conjugate];
%/. zl->1-m/k2+m/k1 //Simplify;
dsigdzl = Simplify[%,{S, α}∈ Reals]
dsigbardz1= 2 dsigdzl /. (S^2)->(1+(cl^2))/4 //Simplify (*Constant multiplied by 2, (S^2) term replaced by (1+cl^2)/2, but
divide by 4 instead to compensate overall 2 *)

Repeat Pion Compton, using polarization sums (2.117 on p 31)

FCClearScalarProducts[];
coef0=-I 4 Pi α; coef1=(k2^2)/(32 Pi (m^2)(k1^2));
{{p1v,p2v,k1v,k2v,k1vr,k2vr}=Momentum/@{p1f,p2f,k1f,k2f,k1r,k2r}};
{ScalarProduct[p1f,p1f]=ScalarProduct[p2f,p2f]=(m^2),
ScalarProduct[k1f,k1f]=ScalarProduct[k2f,k2f]=0,
ScalarProduct[k1r,k1r]=ScalarProduct[k2r,k2r]=0,
ScalarProduct[k1f,k2f]=ScalarProduct[k1r,k2r]=k1 k2 (1-zl),ScalarProduct[p1f,k1f]=ScalarProduct[p1f,k1r]=m k1,ScalarProduct[p1f,k2f]=ScalarProduct[p1f,k2r]=m k2,
ScalarProduct[p1f,p2f]=m (m+k1-k2), ScalarProduct[k1f,p2f]=m k1 - k1 k2 (1-zl), 
ScalarProduct[k2r,p2f]=m k2 -2 k2^2 +k1 k2 (1+zl),
ScalarProduct[k2f,p2f]=m k2 + k1 k2 (1-zl),
ScalarProduct[k1r,p2f]=2 k1^2 + m k1 - k1 k2 (1+zl),
ScalarProduct[k1f,k1r]=2 (k1^2),ScalarProduct[k2f,k2r]=2 (k2^2),
ScalarProduct[k1f,k2r]=ScalarProduct[k1r,k2f]=k1 k2 (1+zl)};

Mmunu=coef0 ((FourVector[k2v,μ]+2 FourVector[p2v,μ])(FourVector[k1v,ν]+2FourVector[p1v,ν])/(2 m k1)-
(2FourVector[p1v,μ]-FourVector[k2v,μ])(2 FourVector[p2v,ν]-FourVector[k1v,ν])/(2 m k2) 
-2 MT[μ,ν])
Mstarrhosogma=ComplexConjugate[Mmunu]/.{μ->ρ,ν->σ}

dsigbardz1=coef1 PolarizationSum[μ,ρ,k2f,k2r]PolarizationSum[ν,σ,k1f,k1r]Mmunu Mstarrhosogma /2 //Contract//Simplify

Electron Compton Scattering, p 55

FCClearScalarProducts[];
{coef0=-I 4 Pi α, coef1=(k2^2)/(32 Pi (m^2)(k1^2))};
{ep1p=Polarization[k1f,Transversality->True],
ep1=Momentum[ep1p],
ep2p=Polarization[k2f,Transversality->True],
ep2=Momentum[ep2p],
ep1ps= ComplexConjugate[ep1p],
ep1s=Momentum[ep1ps],
ep2ps= ComplexConjugate[ep2p],
ep2s=Momentum[ep2ps]};
{p1v,p2v,k1v,k2v}=Momentum/@{p1f,p2f,k1f,k2f};
{{p1sl,p2sl,k1sl,k2sl}=DiracSlash/@{p1f,p2f,k1f,k2f}};
{{ep1sl,ep2ssl}=DiracSlash/@{ep1p,ep2ps}};
{ScalarProduct[p1f,p1f]=ScalarProduct[p2f,p2f]=(m^2),
ScalarProduct[k1f,k1f]=ScalarProduct[k2f,k2f]=0,
ScalarProduct[k1f,k2f]=k1 k2 (1-zl),ScalarProduct[p1f,k1f]=m k1,ScalarProduct[p1f,k2f]=m k2,
ScalarProduct[p1f,p2f]=m (m+k1-k2), ScalarProduct[k1f,p2f]=m k1 - k1 k2 (1-zl), ScalarProduct[k2f,p2f]=m k2 + k1 k2 (1-zl),
ScalarProduct[p1f,ep1p]=ScalarProduct[p1f,ep2p]=0,
ScalarProduct[p1f,ep1ps]=ScalarProduct[p1f,ep2ps]=0,
ScalarProduct[k1f,ep2p]=A,ScalarProduct[k1f,ep2ps]=A\[Conjugate],
ScalarProduct[k2f,ep1p]=B, ScalarProduct[k2f,ep1ps]=B\[Conjugate],
ScalarProduct[p2f,ep1p]=-B, ScalarProduct[p2f,ep1ps]=-B\[Conjugate],
ScalarProduct[p2f,ep2p]=A, ScalarProduct[p2f,ep2ps]=A\[Conjugate],
ScalarProduct[ep1p,ep2ps]=S,ScalarProduct[ep1ps,ep2p]=S\[Conjugate],ScalarProduct[ep1p,ep2p]=T,ScalarProduct[ep1ps,ep2ps]=T\[Conjugate]};

f1 = coef0 SpinorUBar[p2f,m].ep2ssl.(p1sl+k1sl+m).ep1sl.SpinorU[p1f,m]/(2 Pair[p1v,k1v])  //FCI  //DiracSimplify
f2 = coef0 SpinorUBar[p2f,m].ep1sl.(p1sl-k2sl+m).ep2ssl.SpinorU[p1f,m]/(-2 Pair[p1v,k2v])  //FCI //DiracSimplify
f1star=ComplexConjugate[f1] 
f2star=ComplexConjugate[f2]

Sum/average on electron but not photon spin

m11= FermionSpinSum[f1 f1star+ f2 f2star+f1 f2star + f2 f1star]/2; 
%/.DiracTrace->TR //Simplify
%/. zl->1-m/k2+m/k1 //Simplify
Simplify[%,{A,B,S,T,α}∈ Reals]/.T->S
% //Expand
Mbsq=Collect[%,α,Simplify]

Mbsq/. 16 Pi^2 (α^2) -> C //Expand;
Mbsqred=Collect[%,C]/.C-> 16 Pi^2 (α^2)

Differential cross section in lab for massless photon

dsigdzl=coef1 Mbsq;
coef1 Mbsqred

Now sum/average over photon

dsigbardzl= 2 dsigdzl /. (S^2)-> (2-(sl^2))/4 //Simplify ;(*Constant multiplied by 2, (S^2) term replaced by (2-sl^2)/2, but
divide by 4 instead to compensate overall 2 *)
dsigbardzl/( Pi (α^2)(k2^2)/((m^2) (k1^2) )) //Expand;
Collect[C1  %,C1]/. C1->( Pi (α^2)(k2^2)/((m^2) (k1^2) ))

Repeat electron Compton, using polarization sums (2.117 on p 31)

FCClearScalarProducts[];
coef0=-I 4 Pi α; coef1=(k2^2)/(32 Pi (m^2)(k1^2));
{{p1v,p2v,k1v,k2v,k1vr,k2vr}=Momentum/@{p1f,p2f,k1f,k2f,k1r,k2r}};
{{p1sl,p2sl,k1sl,k2sl}=DiracSlash/@{p1f,p2f,k1f,k2f}};
{ScalarProduct[p1f,p1f]=ScalarProduct[p2f,p2f]=(m^2),
ScalarProduct[k1f,k1f]=ScalarProduct[k2f,k2f]=0,
ScalarProduct[k1r,k1r]=ScalarProduct[k2r,k2r]=0,
ScalarProduct[k1f,k2f]=ScalarProduct[k1r,k2r]=k1 k2 (1-zl),ScalarProduct[p1f,k1f]=ScalarProduct[p1f,k1r]=m k1,ScalarProduct[p1f,k2f]=ScalarProduct[p1f,k2r]=m k2,
ScalarProduct[p1f,p2f]=m (m+k1-k2), ScalarProduct[k1f,p2f]=m k1 - k1 k2 (1-zl), 
ScalarProduct[k2r,p2f]=m k2 -2 k2^2 +k1 k2 (1+zl),
ScalarProduct[k2f,p2f]=m k2 + k1 k2 (1-zl),
ScalarProduct[k1r,p2f]=2 k1^2 + m k1 - k1 k2 (1+zl),
ScalarProduct[k1f,k1r]=2 (k1^2),ScalarProduct[k2f,k2r]=2 (k2^2),
ScalarProduct[k1f,k2r]=ScalarProduct[k1r,k2f]=k1 k2 (1+zl)};

Mmunu= coef0 (SpinorUBar[p2f,m].GA[μ].(p1sl+k1sl+m).GA[ν].SpinorU[p1f,m]/(2 Pair[p1v,k1v])+ SpinorUBar[p2f,m].GA[ν].(p1sl-k2sl+m).GA[μ].SpinorU[p1f,m]/(-2 Pair[p1v,k2v]))   //DiracSimplify
Mstarrhosogma=ComplexConjugate[Mmunu]/.{μ->ρ,ν->σ}
 FermionSpinSum[Mmunu Mstarrhosogma/2 //Expand];
Msq=%/.DiracTrace->TR //Simplify;
PolarizationSum[μ,ρ,k2f,k2r]PolarizationSum[ν,σ,k1f,k1r]Msq /2 //Contract//Simplify
%/(32 Pi^2 α^2)
%+1-zl^2//Simplify
%/. zl->1-m/k2+m/k1 //Simplify //Apart
dsigbardz1=coef1 32 Pi^2 α^2 (%+zl^2-1)(*agrees with 2.234 *)
(*PolarizationSum[μ,ρ,k2f,k2r]PolarizationSum[ν,σ,k1f,k1r]Msq /2 //Contract//Simplify /. k2->k1/(1+k1(1-zl)/m) //Simplify
dsigbardz1=  coef1 %*)  (*agrees with 2.234 *)

TOC

QCD Calculations

Color Identities, Table 5.1, p194

SUNT[i]SUNT[i]//SUNSimplify
SUNTrace[SUNT[i,j]]
SUNTrace[SUNT[i,j,k],Explicit->True]
SUNSimplify[SUNTrace[SUNT[i,j,i,l]],Explicit->True]

I SUNSimplify[SUNF[i,j,m]SUNTrace[SUNT[i,j,k]],Explicit->True]

SUNSimplify[SUNF[i,j,a]SUNF[i,j,a]]
SUNSimplify[SUNF[i,j,a]SUNF[i,j,b]]

SUNSimplify[SUNF[i,j,a]SUNF[k,l,a]SUNF[i,k,b]SUNF[j,l,b],Explicit->True]
SUNSimplify[SUNF[i,j,a]SUNF[k,l,a]SUNF[i,j,b]SUNF[k,l,b],Explicit->True]
SUNSimplify[SUNF[i,j,a]SUNF[i,j,a]SUNF[k,l,b]SUNF[k,l,b],Explicit->True]

New Command for version 9, gives more freedom in ordering color flow

t1=SUNTF[{i},a,b]
t1s=ComplexConjugate[t1]
t1 t1s//SUNSimplify

t2=SUNTF[{j},b,c]
t1 t2//SUNSimplify
%  SUNDelta[i,j]
%//SUNDeltaContract
%//SUNSimplify

t1 t2 //SUNSimplify
%  SUNFDelta[a,c]
%//SUNSimplify
t1 t2  SUNFDelta[a,c]
%//SUNSimplify

t3=SUNTF[{i,j},a,b]
%   SUNDelta[i,j]//SUNDeltaContract//SUNSimplify

Reproduce major entries from Table 5.2, p196, sometimes including masses.

(Subscript[q, r])(Subscript[q, s]) -> (Subscript[q, r])(Subscript[q, s]) from Table 5.2, p196

FCClearScalarProducts[];
{ScalarProduct[p1c,p1c]=ScalarProduct[p3c,p3c]=((Subscript[m, r])^2),
ScalarProduct[p2c,p2c]=ScalarProduct[p4c,p4c]=((Subscript[m, s])^2),
ScalarProduct[p1c,p2c]=ScalarProduct[p3c,p4c]=(s-((Subscript[m, r])^2)-((Subscript[m, s])^2))/2,ScalarProduct[p1c,p4c]=ScalarProduct[p2c,p3c]=(-u+((Subscript[m, r])^2)+((Subscript[m, s])^2))/2,
ScalarProduct[p1c,p3c]= ((Subscript[m, r])^2)-t/2,ScalarProduct[p2c,p4c]=((Subscript[m, s])^2)-t/2};

f1=  GlPr[p3c-p1c,μ,i,ν,j]SpinorUBar[p3c,(Subscript[m, r]),3].QGV9[μ,i,a,b].SpinorU[p1c,(Subscript[m, r]),1] SpinorUBar[p4c,(Subscript[m, s]),4].QGV9[ν,j,c,d]. SpinorU[p2c,(Subscript[m, s]),2]//FCI  //Explicit
 f1star=ComplexConjugate[f1]//FCRenameDummyIndices
Mbsq=FermionSpinSum[f1 f1star/qav^2] /.ftrace// Contract  //Simplify 
%/.{(Subscript[m, r])->0} 
%/.{(Subscript[m, s])->0} 
(*TrickMandelstam[%,{s,t,u,2 ((Subscript[m, r])^2)+ 2 ((Subscript[m, s])^2)}] *)

qq -> qq from Table 5.2, p196

FCClearScalarProducts[];
{ScalarProduct[p1d,p1d]=ScalarProduct[p2d,p2d]=(m^2),
ScalarProduct[p3d,p3d]=ScalarProduct[p4d,p4d]=(m^2),
ScalarProduct[p1d,p2d]=ScalarProduct[p3d,p4d]=s/2-(m^2),ScalarProduct[p1d,p4d]=ScalarProduct[p2d,p3d]= (m^2)-u/2,
ScalarProduct[p1d,p3d]=ScalarProduct[p2d,p4d]= (m^2)-t/2};
f1=  GlPr[p3d-p1d,μ,i,ν,j]SpinorUBar[p3d,m].QGV9[μ,i,a,b].SpinorU[p1d,m] SpinorUBar[p4d,m].QGV9[ν,j,c,d]. SpinorU[p2d,m]-
GlPr[p4d-p1d,μ,i,ν,j]SpinorUBar[p4d,m].QGV9[μ,i,c,b].SpinorU[p1d,m] SpinorUBar[p3d,m].QGV9[ν,j,a,d]. SpinorU[p2d,m]//Explicit
 
f1star=ComplexConjugate[f1]//FCRenameDummyIndices
Mbsq=FermionSpinSum[f1 f1star/qav^2] /.ftrace// Contract  //Simplify 
%/. m->0 
%//Apart
9%/4/ (Subscript[g, s])^4//Apart
% (4 (Subscript[g, s])^4/9)

(Subscript[q, r])(Subscript[(OverBar[q]), r]) -> (Subscript[q, s])(Subscript[(OverBar[q]), s]) from Table 5.2, p196

Check mass factor from PDG'08 (39.13)

FCClearScalarProducts[];
{ScalarProduct[p1e,p1e]=ScalarProduct[p2e,p2e]=((Subscript[m, r])^2),
ScalarProduct[p3e,p3e]=ScalarProduct[p4e,p4e]=((Subscript[m, s])^2),
ScalarProduct[p1e,p2e]=(s-2 ((Subscript[m, r])^2))/2,ScalarProduct[p3e,p4e]=(s-2((Subscript[m, s])^2))/2,ScalarProduct[p1e,p4e]=ScalarProduct[p2e,p3e]=(-u+((Subscript[m, r])^2)+((Subscript[m, s])^2))/2,
ScalarProduct[p1e,p3e]=ScalarProduct[p2e,p4e]=(-t+((Subscript[m, r])^2)+((Subscript[m, s])^2))/2};
f1=  GlPr[p2e+p1e,μ,i,ν,j]SpinorVBar[p2e,(Subscript[m, r])].QGV9[μ,i,a,b].SpinorU[p1e,(Subscript[m, r])] SpinorUBar[p3e,(Subscript[m, s])].QGV9[ν,j,c,d]. SpinorV[p4e,(Subscript[m, s])] //Contract //Explicit

f1star=ComplexConjugate[f1]//FCRenameDummyIndices
Mbsq=FermionSpinSum[f1 f1star/qav^2] /.ftrace// Contract  //Simplify 

Mbsq /.{(Subscript[m, r])->0}

GG -> q (OverBar[q]) from Table 5.2, p196

This form of calculation contracts the gluon propagator indices by hand.

Amplitude has two vertices: QGV9 and QGV9 use SMB convention; opposite from built in vertices
For explicit conjugation: First two diagrams ∼(i^3), change sign in conjugate. Third ∼(i^2), doesn't.

FCClearScalarProducts[];
{ScalarProduct[p1g,p1g]=ScalarProduct[p2g,p2g]=0,
ScalarProduct[p3g,p3g]=ScalarProduct[p4g,p4g]=(m^2),
ScalarProduct[p1g,p2g]=s/2,
ScalarProduct[p3g,p4g]=(s-2(m^2))/2,ScalarProduct[p1g,p4g]=ScalarProduct[p2g,p3g]=(-u+(m^2))/2,
ScalarProduct[p1g,p3g]=ScalarProduct[p2g,p4g]=(-t+(m^2))/2};
f1= SpinorUBar[p3g,m].QGV9[μ,i,a,b].QPr[p3g-p1g,m].QGV9[ν,j,b,c].SpinorV[p4g,m]//FCI//Explicit
(*f1star=ComplexConjugate[f1]//FCRenameDummyIndices;*)
f2= SpinorUBar[p3g,m].QGV9[ν,j,a,b].QPr[p3g-p2g,m].QGV9[μ,i,b,c].SpinorV[p4g,m]//FCI//Explicit
(*f2star=ComplexConjugate[f2]//FCRenameDummyIndices;*)
f3= SpinorUBar[p3g,m].QGV9[α,k,a,c].SpinorV[p4g,m](-ScPr[p1g+p2g])GV9[{p1g,μ,i},{p2g,ν,j},{-p1g-p2g,α,k}]//FCI//Explicit
(*f3star=ComplexConjugate[f3]//FCRenameDummyIndices*)
ft=f1+f2+f3;
ftstar=ComplexConjugate[ft]/.{μ->ρ,ν->σ}//FCRenameDummyIndices;
(*ftstar2=f1star+f2star+f3star;*)
POLs= PolarizationSum[μ,ρ,p1g,p2g]PolarizationSum[ν,σ,p2g,p1g];
(*Define polarization sums using p1r=p2 and p2r=p1 in CM*)

gval=FermionSpinSum[ft ftstar  /Gav^2//Expand];
eval=(POLs gval//Contract)/.ftrace; 
eval1=eval/(Subscript[g, s])^4//Simplify
eval5=TrickMandelstam[eval1,{s,t,u,2 m^2}]
eval5/.{s->4 m^2, t->-m^2, u-> -m^2}  (*PDG: 7/48*)
eval2=eval1/. m->0
eval3=TrickMandelstam[eval2,{s,t,u,0}];
spole= 3(t^2+u^2)/(8 s^2);
eval4=eval3+spole//Combine//Simplify;
(Subscript[g, s])^4(TrickMandelstam[eval4,{s,t,u,0}]-spole)
(*eval5//InputForm*)

Repeat, using full Feynman gauge gluon propagator

f3a= SpinorUBar[p3g,m].QGV9[α1,k1,a,c].SpinorV[p4g,m](GlPr[p1g+p2g,α1,k1,α2,k2])GV9[{p1g,μ,i},{p2g,ν,j},{-p1g-p2g,α2,k2}]//FCI//Explicit
(*f3star=ComplexConjugate[f3]//FCRenameDummyIndices*)
fta=f1+f2+f3a;
ftastar=ComplexConjugate[fta]/.{μ->ρ,ν->σ}//FCRenameDummyIndices;
(*ftstar2=f1star+f2star+f3star;*)

gval=FermionSpinSum[fta ftastar  /Gav^2//Expand];
eval=(POLs gval//Contract)/.ftrace; 
eval1=eval/(Subscript[g, s])^4//Simplify
eval5=TrickMandelstam[eval1,{s,t,u,2 m^2}]
eval5/.{s->4 m^2, t->-m^2, u-> -m^2}  (*PDG: 7/48*)
eval2=eval1/. m->0
eval3=TrickMandelstam[eval2,{s,t,u,0}];
spole= 3(t^2+u^2)/(8 s^2);
eval4=eval3+spole//Combine//Simplify;
(Subscript[g, s])^4(TrickMandelstam[eval4,{s,t,u,0}]-spole)

GG -> GG  from Table 5.2, p196

First method: contract unnecessary indices.

FCClearScalarProducts[];
{ScalarProduct[p1h,p1h]=ScalarProduct[p2h,p2h]=0,
ScalarProduct[p3h,p3h]=ScalarProduct[p4h,p4h]=0,
ScalarProduct[p1h,p2h]=ScalarProduct[p3h,p4h]=s/2,
ScalarProduct[p1h,p4h]=ScalarProduct[p2h,p3h]=(-u)/2,
ScalarProduct[p1h,p3h]=ScalarProduct[p2h,p4h]=(-t)/2};
(*prodproc[i_,j_]:=flist[[i]] fliststar[[j]]//Explicit
prodcol[i_,j_]:=CalcColorFactor[prodproc[i,j]]
prodcont[i_,j_]:=POLs prodcol[i,j]/. D->4 //Contract //Simplify*)

f1=GV49a[{p1h,μ,i},{p2h,ν,j},{-p3h,ρ,l},{-p4h,σ,k}]//Explicit
f2=GV9[{p1h,μ,i},{p2h,ν,j},{-p1h-p2h,α,a}](-ScPr[p1h+p2h])GV9[{-p3h,ρ,l},{-p4h,σ,k},{p3h+p4h,α,a}]//Explicit
f3=GV9[{p1h,μ,i},{-p3h,ρ,l},{-p1h+p3h,α,a}](-ScPr[p1h-p3h])GV9[{p2h,ν,j},{-p4h,σ,k},{-p2h+p4h,α,a}]//Explicit
f4=GV9[{p1h,μ,i},{-p4h,σ,k},{-p1h+p4h,α,a}](-ScPr[p1h-p4h])GV9[{p2h,ν,j},{-p3h,ρ,l},{-p2h+p3h,α,a}]//Explicit
ft=f1+f2+f3+f4;

ftstar=ComplexConjugate[ft]/.{μ->μp,ν->νp,ρ->ρp,σ->σp}//FCRenameDummyIndices;
POLs= PolarizationSum[μ,μp,p1h,p2h]PolarizationSum[ν,νp,p2h,p1h] PolarizationSum[ρ,ρp,p3h,p4h]PolarizationSum[σ,σp,p4h,p3h];

prod= ft ftstar //CalcColorFactor;
rprod=POLs prod//Contract;
Mbsqf=rprod/Gav^2 //Simplify

redampsq=TrickMandelstam[Mbsqf,{s,t,u,0}]4/(9(Gstrong^4))//Simplify
(*verify that results agrees with Table 5.2*)altform=(s^2+u^2)/t^2+(s^2+t^2)/u^2+(t^2+u^2)/s^2+3
redampsq-altform 
TrickMandelstam[%,{s,t,u,0}]

Second method: use full gluon propagators (runs almost as fast)

FCClearScalarProducts[];
{ScalarProduct[p1h,p1h]=ScalarProduct[p2h,p2h]=0,
ScalarProduct[p3h,p3h]=ScalarProduct[p4h,p4h]=0,
ScalarProduct[p1h,p2h]=ScalarProduct[p3h,p4h]=s/2,
ScalarProduct[p1h,p4h]=ScalarProduct[p2h,p3h]=(-u)/2,
ScalarProduct[p1h,p3h]=ScalarProduct[p2h,p4h]=(-t)/2};
(*prodproc[i_,j_]:=flist[[i]] fliststar[[j]]//Explicit
prodcol[i_,j_]:=CalcColorFactor[prodproc[i,j]]
prodcont[i_,j_]:=POLs prodcol[i,j]/. D->4 //Contract //Simplify*)
f1=GV49a[{p1h,μ,i},{p2h,ν,j},{-p3h,ρ,l},{-p4h,σ,k}]//Explicit
f2=GV9[{p1h,μ,i},{p2h,ν,j},{-p1h-p2h,α1,a1}]GlPr[p1h+p2h,α1,a1,α2,a2]GV9[{-p3h,ρ,l},{-p4h,σ,k},{p3h+p4h,α2,a2}]//Explicit
f3=GV9[{p1h,μ,i},{-p3h,ρ,l},{-p1h+p3h,α1,a1}]GlPr[p1h-p3h,α1,a1,α2,a2]GV9[{p2h,ν,j},{-p4h,σ,k},{-p2h+p4h,α2,a2}]//Explicit
f4=GV9[{p1h,μ,i},{-p4h,σ,k},{-p1h+p4h,α1,a1}]GlPr[p1h-p4h,α1,a1,α2,a2]GV9[{p2h,ν,j},{-p3h,ρ,l},{-p2h+p3h,α2,a2}]//Explicit
ft=f1+f2+f3+f4;

ftstar=ComplexConjugate[ft]/.{μ->μp,ν->νp,ρ->ρp,σ->σp}//FCRenameDummyIndices;
POLs= PolarizationSum[μ,μp,p1h,p2h]PolarizationSum[ν,νp,p2h,p1h] PolarizationSum[ρ,ρp,p3h,p4h]PolarizationSum[σ,σp,p4h,p3h];
prod= ft ftstar //CalcColorFactor;
rprod=POLs prod//Contract;
Mbsqf=rprod/Gav^2 //Simplify

redampsq=TrickMandelstam[Mbsqf,{s,t,u,0}]4/(9(Gstrong^4))//Simplify
(*verify that results agrees with Table 5.2*)altform=(s^2+u^2)/t^2+(s^2+t^2)/u^2+(t^2+u^2)/s^2+3
redampsq-altform 
TrickMandelstam[%,{s,t,u,0}]

Third method: work with helicity amplitudes

Use circular polarization vectors defined on p31. Reflection and crossing relate many amplitudes. In following, we will ignore overall sign of each amplitude. In particular, there are only 6 nonvanishing helicity amplitudes. This follows from Parke-Taylor theorem. In convention in which all 4 are outgoing or incoming
the amplitudes with 4 same or 3 same+1 different helicities vanish. e.g. RRLL->rrrr in convention of all
incoming. This vanishes, as does RRRL->rrlr and RLRR-> rlll.

 3 are related by reflection to the other 3. The three remaining are related by crossing symmetry.

let li=+1 (right) or -1 (left)

prodproc[i_,j_]:=flist[[i]] fliststar[[j]]//Explicit
prodcol[i_,j_]:=CalcColorFactor[prodproc[i,j]]
rodcont[i_,j_]:=POLs prodcol[i,j]/. D->4 //Contract //Simplify

FCClearScalarProducts[];
{ep1p=Polarization[p1h,Transversality->True],
ep2p=Polarization[p2h,Transversality->True],
ep3p=Polarization[p3h,Transversality->True],
ep4p=Polarization[p4h,Transversality->True]};  
(*ep3p=Polarization[p3h,Transversality->True] //ComplexConjugate,
ep4p=Polarization[p4h,Transversality->True]//ComplexConjugate}*)
{ScalarProduct[p1h,p1h]=ScalarProduct[p2h,p2h]=0,
ScalarProduct[p3h,p3h]=ScalarProduct[p4h,p4h]=0,
ScalarProduct[p1h,p2h]=ScalarProduct[p3h,p4h]=s/2,
ScalarProduct[p1h,p4h]=ScalarProduct[p2h,p3h]=(-u)/2,
ScalarProduct[p1h,p3h]=ScalarProduct[p2h,p4h]=(-t)/2};
(*{cz=(t-u)/s, sz=2Sqrt[t u]/s,
ScalarProduct[ep1p,ep2p]=-(1+l1 l2)/2,
ScalarProduct[ep1p,ep3p]=-(cz+l1 l3)/2,
ScalarProduct[ep1p,ep4p]=-(cz-l1 l4)/2,
ScalarProduct[ep2p,ep3p]=-(cz-l2 l3)/2,
ScalarProduct[ep2p,ep4p]=-(cz+l2 l4)/2,
ScalarProduct[ep3p,ep4p]=-(1+l3 l4)/2};*)
{cz=(t-u)/s, sz=2Sqrt[t u]/s,
ScalarProduct[ep1p,ep2p]=+(1+l1 l2)/2,
ScalarProduct[ep1p,ep3p]=-(cz+l1 l3)/2,  
ScalarProduct[ep1p,ep4p]=+(cz-l1 l4)/2,
ScalarProduct[ep2p,ep3p]=+(cz-l2 l3)/2,
ScalarProduct[ep2p,ep4p]=-(cz+l2 l4)/2,
ScalarProduct[ep3p,ep4p]=+(1+l3 l4)/2};
(*old version had opposite signs for eps1,3; this form consistent with SMB*)
(*formulae are really for conjugates of 3 and 4*)

(*{ScalarProduct[ep1p,p2h]=ScalarProduct[ep2p,p1h]=ScalarProduct[ep3p,p4h]=ScalarProduct[ep4p,p3h]=0,
ScalarProduct[ep1p,p3h]=ScalarProduct[ep2p,p3h]=ScalarProduct[ep3p,p2h]=ScalarProduct[ep4p,p2h]=-Sqrt[s/2]sz/2,
ScalarProduct[ep1p,p4h]=ScalarProduct[ep2p,p4h]=ScalarProduct[ep3p,p1h]=ScalarProduct[ep4p,p1h]=Sqrt[s/2]sz/2};*)
{ScalarProduct[ep1p,p2h]=ScalarProduct[ep2p,p1h]=ScalarProduct[ep3p,p4h]=ScalarProduct[ep4p,p3h]=0,
ScalarProduct[ep1p,p3h]=ScalarProduct[ep2p,p4h]=ScalarProduct[ep3p,p2h]=ScalarProduct[ep4p,p1h]=Sqrt[s/2]sz/2,
ScalarProduct[ep1p,p4h]=ScalarProduct[ep2p,p3h]=ScalarProduct[ep3p,p1h]=ScalarProduct[ep4p,p2h]=-Sqrt[s/2]sz/2};

f1=GV49a[{p1h,μ,i},{p2h,ν,j},{-p3h,ρ,l},{-p4h,σ,k}]//Explicit
f2=GV9[{p1h,μ,i},{p2h,ν,j},{-p1h-p2h,α,a}]GlPr[p1h+p2h,α,a,β,b]GV9[{-p3h,ρ,l},{-p4h,σ,k},{p3h+p4h,β,b}]//Explicit
f3=GV9[{p1h,μ,i},{-p3h,ρ,l},{-p1h+p3h,α,a}]GlPr[p1h-p3h,α,a,β,b]GV9[{p2h,ν,j},{-p4h,σ,k},{-p2h+p4h,β,b}]//Explicit
f4=GV9[{p1h,μ,i},{-p4h,σ,k},{-p1h+p4h,α,a}]GlPr[p1h-p4h,α,a,β,b]GV9[{p2h,ν,j},{-p3h,ρ,l},{-p2h+p3h,β,b}]//Explicit
flist=List[f1,f2,f3,f4]

(*POLp=Pair[Momentum[ep1p],LorentzIndex[μ]]Pair[Momentum[ep2p],LorentzIndex[ν]]ComplexConjugate[Pair[Momentum[ep3p],LorentzIndex[ρ]]]ComplexConjugate[Pair[Momentum[ep4p],LorentzIndex[σ]]]
*)
POLp=Pair[Momentum[ep1p],LorentzIndex[μ]]Pair[Momentum[ep2p],LorentzIndex[ν]]Pair[Momentum[ep3p],LorentzIndex[ρ]]Pair[Momentum[ep4p],LorentzIndex[σ]]
famplist={a1,a2,a3,a4};
Do[famplist[[iv]]=POLp flist[[iv]]//Contract//Simplify,{iv,4}];
famplist;  
famp=Sum[famplist[[iv]],{iv,1,4}]//SUNSimplify//Simplify

Print["-------------------------------------"]
Print["-----------RRRL (should vanish)--------"]

famp /.{l1->1,l2->1,l3->1,l4->-1};     (*verify RRRL vanishes*)
SUNSimplify[%,SUNFJacobi->True]//Simplify
TrickMandelstam[%,{s,t,u,0}]

Print["-------------------------------------"]
Print["-----------RRLL (should vanish)--------"]

famp /.{l1->1,l2->1,l3->-1,l4->-1};   (*verify RRLL vanishes*)
SUNSimplify[%,SUNFJacobi->True]//Simplify
TrickMandelstam[%,{s,t,u,0}]

Print["-------------------------------------"]
Print["-----------RRRR--------"]

famp /.{l1->1,l2->1,l3->1,l4->1};    (* RRRR *)
SUNSimplify[%,SUNFJacobi->True]//Simplify;
fRRRR=TrickMandelstam[%,{s,t,u,0}] (*shows expect symmetry under t <-> u and k <-> l*)
fRRRR-(fRRRR/.{t:>u,u:>t, k:>l, l:>k})(*this shows the symmetry explicitly*)
(*hand calculation shows that fRRRR=-2 i g^2 (bik bjl s/u-bil bjk s/t)*)

Print["-------------------------------------"]
Print["-----------LLLL--------"]

famp /.{l1->-1,l2->-1,l3->-1,l4->-1};  (* LLLL *)
SUNSimplify[%,SUNFJacobi->True]//Simplify;
fLLLL=TrickMandelstam[%,{s,t,u,0}] (*shows expected reflection invariance*)
fRRRR-fLLLL

Print["-------------------------------------"]
Print["-----------RLRL--------"]

famp /.{l1->1,l2->-1,l3->1,l4->-1};  (* RLRL *)
SUNSimplify[%,SUNFJacobi->True]//Simplify;
fRLRL=TrickMandelstam[%,{s,t,u,0}]

Print["-------------------------------------"]
Print["-----------LRLR--------"]

famp /.{l1->-1,l2->1,l3->-1,l4->+1};(* LRLR *)
SUNSimplify[%,SUNFJacobi->True]//Simplify;
fLRLR=TrickMandelstam[%,{s,t,u,0}]

Print["-------------------------------------"]
Print["-----------RLLR--------"]

famp /.{l1->1,l2->-1,l3->-1,l4->+1};(* RLLR *)
SUNSimplify[%,SUNFJacobi->True]//Simplify;
fRLLR=TrickMandelstam[%,{s,t,u,0}] (*shows t-u and k-l relation to fLRLR*)

flist=List[fRRRR,fRLRL,fRLLR]
fliststar={fRRRRs,fRLRLs,fRLLRs}=ComplexConjugate/@{fRRRR,fRLRL,fRLLR}//FCRenameDummyIndices

CalcColorFactor[2 fRRRR fRRRRs/(Gav^2) //Explicit];
(*Multiply by 2 for opposite parity states.Results have expected s-t-u crossing symmetry*)
msqRRRR=TrickMandelstam[%,{s,t,u,0}] 
msqRLLR =CalcColorFactor[2fRLLR fRLLRs /(Gav^2)//Explicit]//Simplify
msqRLRL=CalcColorFactor[2fRLRL fRLRLs /(Gav^2)//Explicit]//Simplify
Mbsqf=msqRRRR+msqRLLR +msqRLRL//Simplify
redampsq=TrickMandelstam[Mbsqf,{s,t,u,0}]4/(9(Gstrong^4))//Simplify
(*verify that results agrees with Table 5.2*)
(*%//Apart
TrickMandelstam[%,{s,t,u,0}]*)
altform=(s^2+u^2)/t^2+(s^2+t^2)/u^2+(t^2+u^2)/s^2+3
redampsq-altform 
TrickMandelstam[%,{s,t,u,0}]

q (OverBar[q]) →  (Subscript[q, 0]) (Subscript[(OverBar[q]), 0]) including scalar mass (QCD only, no gluino)

SetMandelstam[s, t, u, p1qs, p2qs,- p3qs, -p4qs,0,0,m,m];
fqs=Spinor[-p2qs].QGV9[μ,i,a,b].Spinor[p1qs]  GlPr[p1qs+p2qs,μ,i,ν,j] QsGV9[p3qs-p4qs,ν,j,c,d];
(*fqsstar=-ComplexConjugate[fqs]  /.{μ-> ρ, i->k,ν-> σ, j->l}*) (*explicit sign since I in QGV not conjugated*)
fqsstar=ComplexConjugate[fqs]//FCRenameDummyIndices;
Mbsq=FermionSpinSum[fqs fqsstar/qav^2]/.ftrace//Simplify
TrickMandelstam[%,{s,t,u,2(m^2)}]

GG→  (Subscript[q, 0]) (Subscript[(OverBar[q]), 0]) (massless scalar)

SetMandelstam[s, t, u, p1gs, p2gs,- p3gs, -p4gs,0,0,0,0];
f1=GV9[{p1gs,μ,i},{p2gs,ν,j},{-p1gs-p2gs,α,m}]GlPr[p1gs+p2gs,α,m,β,n]QsGV9[p3gs-p4gs,β,n,a,b]//Contract//Simplify//SUNDeltaContract;
f2=QsQsGGV9[{μ,i},{ν,j},a,b];
f3=(QsGV9[2 p3gs-p1gs,μ,i,a,c]).(QsGV9[ p3gs-p1gs-p4gs,ν,j,c,b])ScPr[p3gs-p1gs];
f4=(QsGV9[ p3gs+p1gs-p4gs,ν,j,a,c]).(QsGV9[-2 p4gs+p1gs,μ,i,c,b])ScPr[p4gs-p1gs];
ftot=f1+f2+f3+f4;
ftotstar=ComplexConjugate[ftot]/.{μ->ρ,ν->σ}//FCRenameDummyIndices;
POLs= PolarizationSum[μ,ρ,p1gs,p2gs]PolarizationSum[ν,σ,p2gs,p1gs];
m1sq=ftot.ftotstar//CalcColorFactor;
m2sq=POLs m1sq/Gav^2//Contract//Simplify
m3sq=TrickMandelstam[m2sq,{s,t,u,0}]
%-Gstrong^4 (3 ((u-t)/s)^2/32+7/96)//Simplify  (*verify*)
TrickMandelstam[%,{s,t,u,0}]

TOC

Fermi Theory

μ Decay, p 250

FCClearScalarProducts;
Clear[m]
masslesse ={(Subscript[p, e])-> (Subscript[Ε, e]), (Subscript[m, e])->0};
{coef0=-4 I (Subscript[G, F])/Sqrt[2],coef1= (Subscript[p, e])/(m 2^8 Pi^4)};
(*{(Subscript[p, e])=Sqrt[(Subscript[Ε, e]) -((Subscript[m, e])^2)]};   don't enforce yet*)
{{p1sl,p2sl,p3sl,p4sl, ssl,sesl}=DiracSlash/@{p1j,p2j,p3j,p4j,sm,se}};
{p1v,p2v,p3v,p4v,sv,sev}=Momentum/@{p1j,p2j,p3j,p4j,sm,se};
{ScalarProduct[sm,p1j]=0, ScalarProduct[p1j,p1j]=(m^2),
ScalarProduct[p1j,p4j]=m (Subscript[Ε, e]),
ScalarProduct[p2j,p2j]=ScalarProduct[p3j,p3j]=0,
ScalarProduct[p4j,p4j]=((Subscript[m, e])^2),ScalarProduct[sm,sm]=ScalarProduct[se,se]=-1,ScalarProduct[se,p4j]=0,
ScalarProduct[sm,p4j]=-(Subscript[p, e]) (Subscript[c, θ]),
ScalarProduct[se,p1j]=m (Subscript[p, e])(Subscript[c, he])/(Subscript[m, e]),
ScalarProduct[se,sm]= (1-(Subscript[Ε, e])/(Subscript[m, e]))  (Subscript[c, θ]) (Subscript[c, he]) - (Subscript[c, ss])}(* (Subscript[c, θ]) = cos of angle between pe and muon spin.
(Subscript[c, he]) = cos of angle between pe and electron spin (2 x helicity).
 (Subscript[c, ss])= cos of angle between muon and electron spin vectors in their respective rest frames.
Expressions for se follow from (2.176) and ensure smooth (Subscript[m, e])->0 limit *)

First, use Mbarsq calculated by hand

t13=2 GA[μ].p1sl.GA[ν].p3sl.GA[6]
t24=2 GA[μ].p2sl.GA[ν].p4sl.GA[6]
Mbarsq= ((Subscript[G, F])^2)TR[t13] TR[t24]//Contract//Simplify (*Define Mbarsq by hand*)

Now, allow muon spin polarization, by hand (SpinPolarizationSum option resets all polarization sums, not just the muon's)

Pmus=(1+ GA[5].ssl)/2
(Subscript[M, fi])=coef0 SpinorUBar[p3j,0].GA[μ].GA[7].Pmus.SpinorU[p1j,m] SpinorUBar[p4j,(Subscript[m, e])].GA[μ].GA[7]. SpinorV[p2j,0] //FCI
(Subscript[Mstar, fi])=ComplexConjugate[(Subscript[M, fi])] //FCRenameDummyIndices
(Subscript[Mbarsq, fi]) = (FermionSpinSum[(Subscript[M, fi]) (Subscript[Mstar, fi])]) 
%/.DiracTrace->TR// Contract //Simplify
Mbarsq= % /.masslesse //Simplify

Recalculate in form to allow integration over neutrino momenta. Also, include electron polarization.

Pmus=(1+ GA[5].ssl)/2;
Pes=(1+ GA[5].sesl)/2;
(Subscript[M, fi])=coef0 SpinorUBar[p3j,0].GA[μ].GA[7].Pmus.SpinorU[p1j,m] SpinorUBar[p4j,(Subscript[m, e])].Pes.GA[μ].GA[7]. SpinorV[p2j,0] //FCI
(Subscript[Mstar, fi])=ComplexConjugate[(Subscript[M, fi])] //FCRenameDummyIndices
(Subscript[Mbarsq, fi]) = (FermionSpinSum[(Subscript[M, fi]) (Subscript[Mstar, fi])]) /. {p2sl->GA[α], p3sl-> GA[β]}
Mfull= %/.DiracTrace->TR// Contract //Simplify
(*Mbarsq= % /.masslesse //Simplify*) (*me->0 limit is singular; must treat separately*)
Mfull + (Mfull/.se->-se)//Simplify

qv=p1v-p4v;
(Subscript[T, αβ])=Pi(2 FV[qv,α]FV[qv,β]+Pair[qv,qv]MT[α,β])/6//ExpandScalarProduct
(* Use (6.46) *)
Intfull=Contract[(Subscript[T, αβ]) Mfull] //Simplify

Intval=Simplify[% /.masslesse ,  (Subscript[Ε, e])>0]/.  (Subscript[Ε, e])-> ϵ m/2;
Collect[%,m]

Massless case

dgamvaldedcos= coef1 (m/2) Intval/. (Subscript[p, e])-> ϵ m/2;
Collect[%,m]
%/. (Subscript[c, he]) -> -1  (*simulates sum of  (Subscript[c, he])-1 over (Subscript[c, he])= +1 and -1*)
Integrate[%,{(Subscript[c, θ]),-1,1}]
Integrate[%,{ϵ,0,1}]

Massive case

dgamfulldedcos= coef1 Intfull
(% /. {(Subscript[c, he])->+1, (Subscript[c, ss])->+1}) + (% /. {(Subscript[c, he])->-1,  (Subscript[c, ss])->-1})//Simplify;
((% /. (Subscript[c, θ])->+1) + (% /. (Subscript[c, θ])->-1))/2//Simplify
 m %/2 /.{(Subscript[p, e])-> ϵ m/2, (Subscript[Ε, e]) -> ϵ m/2, (Subscript[m, e])-> 0}
Integrate[%,{(Subscript[c, θ]),-1,1}]
Integrate[%,{ϵ,0,1}]

dred=dgamfulldedcos/(((Subscript[G, F])^2)(Subscript[p, e]) /(48 Pi^3 m))//Simplify
%/.masslesse//Simplify  (*check massless limit*)
(*%/.(Subscript[c, he])->+1;
Coefficient[dred,(Subscript[c, ss])]
Coefficient[dred,(Subscript[c, he])]
Coefficient[%, (Subscript[c, θ])]
dred/.{(Subscript[c, ss])->0, (Subscript[c, θ])->0,(Subscript[c, he])->0}*)   (*remove angular and spin terms*)

Full form agrees with Commins and Bucksbaum (3.29)

β Decay (SuperMinus[e]) Angular Distribution, 6.124, using trace techniques

Calculate angular distribution only, not absolute rate. Neglect proton momentum. 1=n, 2=p, 3=e, 4=ν.

FCClearScalarProducts;
{ScalarProduct[sn,p1k]=0, ScalarProduct[sn,sn]=-1,ScalarProduct[sn,p3k]=-(Subscript[β, e]) (Subscript[e, e])(Subscript[c, en]),ScalarProduct[sn,p4k]=-(Subscript[e, ν]) (Subscript[c, νn]),
ScalarProduct[p1k,p1k]=ScalarProduct[p2k,p2k]=ScalarProduct[p1k,p2k]=(m^2),
ScalarProduct[p1k,p3k]=ScalarProduct[p2k,p3k]=m (Subscript[e, e]),
ScalarProduct[p1k,p4k]=ScalarProduct[p2k,p4k]=m (Subscript[e, ν]),
ScalarProduct[p3k,p3k]=((Subscript[m, e])^2),ScalarProduct[p4k,p4k]=0,ScalarProduct[p3k,p4k]=(Subscript[e, ν]) (Subscript[e, e]) (1-(Subscript[β, e])(Subscript[c, eν]))};

Pns=(1+ GA[5].GS[sn])/2;
f1=(Subscript[G, F]) Spinor[p2k,m].GA[μ].(1-λ GA[5]).Pns.Spinor[p1k,m] Spinor[p3k,(Subscript[m, e])].GA[μ].GA[7].Spinor[-p4k,0]/(2 Sqrt[2])
f1star=ComplexConjugate[f1]/. μ-> ν  (*could use //FCRenameDummyIndices*)
Msq= FermionSpinSum[f1 f1star]/.DiracTrace->TR// Contract //Simplify
Msq1=%/.{(Subscript[c, eν])->0,(Subscript[c, νn])->0 }//Simplify  (*average over ν direction*)

TOC

Standard EW Model

Higgs decays from p369

Decay to two  f (OverBar[f])

FCClearScalarProducts;
{Coef=(Subscript[C, f])  (Subscript[p, f]) /(8 Pi ((Subscript[M, H])^2)),ScalarProduct[p1l]=ScalarProduct[p2l]=(m^2),ScalarProduct[p1l,p2l]=((Subscript[M, H])^2)/2-(m^2)};

M1=-I m Sqrt[Sqrt[2] (Subscript[G, F])]Spinor[p1l,m].Spinor[-p2l,m]
M1star=ComplexConjugate[M1]
Mbsq= FermionSpinSum[M1 M1star//Explicit]/. DiracTrace->TR 
Gam=Coef Mbsq
% /((m^2)  (Subscript[M, H]))/.{ (Subscript[p, f]) -> (Subscript[β, f]) (Subscript[M, H])/2, m->(Subscript[M, H])Sqrt[1- ((Subscript[β, f])^2)]/2}//Simplify
% ((m^2)  (Subscript[M, H]))

Decay to (SuperPlus[W])(SuperMinus[W])

FCClearScalarProducts;
{Coef=(Subscript[C, f])   (Subscript[p, f]) /(8 Pi ((Subscript[M, H])^2)),(Subscript[C, f]) =1, ScalarProduct[p1m,p1m]=ScalarProduct[p2m,p2m]=(m^2),ScalarProduct[p1m,p2m]=((Subscript[M, H])^2)/2-(m^2)};
{p1v,p2v}=Momentum/@{p1m,p2m};

Mmunu=2 I MT[μ,ν] Sqrt[Sqrt[2] (Subscript[G, F])] (m^2)
Mmunustar=-2 I MT[ρ,σ] Sqrt[Sqrt[2] (Subscript[G, F])] (m^2)
Msq= Mmunu Mmunustar
(*Polsum1=(-MT[μ,ρ] +FourVector[p1m,μ]FourVector[p1m,ρ]/(m^2))
Polsum2=(-MT[ν,σ] +FourVector[p2m,ν]FourVector[p2m,σ]/(m^2))
Msqfin=Polsum1 Polsum2 Msq  //Contract  //Simplify*)
(* Above code does polarization sums by hand. Below uses built in routines. *)
Msqfin=PolarizationSum[μ,ρ,p1m] PolarizationSum[ν,σ,p2m]   Msq  //Contract//Simplify
% /. m->Sqrt[x] (Subscript[M, H])/2 //Simplify
Gamm=Coef % /. (Subscript[p, f])->Sqrt[(1-x)](Subscript[M, H])/2

(SuperMinus[e])(SuperPlus[e])  → (SuperMinus[W])(SuperPlus[W]), figure 7.30, p360

Ignore fermion masses. Calculate (SuperMinus[e])(L) (SuperPlus[e]) (R) and RL separately.

FCClearScalarProducts;
SetMandelstam[s, t, u, p1ee, p2ee,- p3ee, -p4ee,0,0, (Subscript[M, W]),(Subscript[M, W])];
{pi=Sqrt[s]/2, pf=Sqrt[s- 4 ((Subscript[M, W])^2)]/2,Ef=Sqrt[s]/2,coefww=pf/(32 Pi s pi)};
ksub={t-> ((Subscript[M, W])^2)-2 pi (Ef-pf (Subscript[c, θ])),u-> ((Subscript[M, W])^2)-2 pi (Ef+pf (Subscript[c, θ]))};
gsub={g->2  Sqrt[Sqrt[2] (Subscript[G, F])]  (Subscript[M, W])}; 
(*numsub={ (Subscript[G, F])->1.17 10^(-5), (Subscript[M, W])->80.4, (Subscript[M, Z])-> 91.19, hbcsq->(.197)^2 10^(10) (*pb-GeV^2*), (Subscript[c, w])-> 0.88,(Subscript[s, w])->0.48}*)
numsub={ (Subscript[G, F])->1.17 10^(-5), (Subscript[M, W])->80.4, (Subscript[M, Z])-> 91.19, hbcsq->(.197)^2 10^(10) (*pb-GeV^2*), cw-> 0.88,sw->0.48}

Include all three diagrams

fnuL=Spinor[-p2ee].wfbfL[ν].QPr[p1ee-p3ee,0].wfbfL[μ].Spinor[p1ee]
fAL=Spinor[-p2ee].afbfL[α,-1].Spinor[p1ee]EwPr[p1ee+p2ee,0]aww[-p3ee, μ, p1ee+p2ee,α, -p4ee, ν]
fZL=Spinor[-p2ee].zfbfL[α,-1/2,-1].Spinor[p1ee]EwPr[p1ee+p2ee,(Subscript[M, Z])]zww[-p3ee, μ, p1ee+p2ee,α, -p4ee, ν]
flistL=List[fnuL,fAL,fZL]
fliststarL={fnuLs,fALs,fZLs}=ComplexConjugate/@{fnuL,fAL,fZL}/.{μ-> ρ, ν-> σ,α->β}//Explicit
fnuR=0
fAR=Spinor[-p2ee].afbfR[α,-1].Spinor[p1ee]EwPr[p1ee+p2ee,0]aww[-p3ee, μ, p1ee+p2ee,α,-p4ee, ν]
fZR=Spinor[-p2ee].zfbfR[α,-1/2,-1].Spinor[p1ee]EwPr[p1ee+p2ee,(Subscript[M, Z])]zww[-p3ee, μ, p1ee+p2ee,α, -p4ee, ν]
flistR=List[fnuR,fAR,fZR]
fliststarR={fnuRs,fARs,fZRs}=ComplexConjugate/@{fnuR,fAR,fZR}/.{μ-> ρ, ν-> σ,α->β}//Explicit

Lsqpol=Table[0,{iv,1,3},{jv,1,3}];
Rsqpol=Table[0,{iv,1,3},{jv,1,3}];
Do[Lsqpol[[iv,jv]] =((FermionSpinSum[PolarizationSum[μ,ρ,p3ee] PolarizationSum[ν,σ,p4ee] flistL[[iv]] fliststarL[[jv]]/eav^2//Expand//EpsEvaluate //Explicit]//DiracTrick)
/.DiracTrace->TR  //Contract)/.LeviCivita[][p1ee,p2ee,p3ee,p4ee]->0   //Simplify;
Print[iv,jv,":  ",Lsqpol[[iv,jv]]],{iv,1,3},{jv,1,3}]
Do[Rsqpol[[iv,jv]] =((FermionSpinSum[PolarizationSum[μ,ρ,p3ee] PolarizationSum[ν,σ,p4ee] flistR[[iv]] fliststarR[[jv]]/eav^2//Expand//EpsEvaluate //Explicit]//DiracTrick)
/.DiracTrace->TR  //Contract)/.LeviCivita[][p1ee,p2ee,p3ee,p4ee]->0  //Simplify;
Print[iv,jv,":  ",Rsqpol[[iv,jv]]],{iv,1,3},{jv,1,3}]

Mnusqpol=Lsqpol[[1,1]]+Rsqpol[[1,1]] //Simplify
MnuAsqpol=Sum[Lsqpol[[iv,jv]]+Rsqpol[[iv,jv]] //Simplify,{iv,1,2},{jv,1,2}]//Simplify
(*MnuAZsqpol=Sum[Lsqpol[[iv,jv]] +Rsqpol[[iv,jv]]//Simplify,{iv,1,3},{jv,1,3}]//Simplify *)
MnuAZsqpolL=Sum[Lsqpol[[iv,jv]] //Simplify,{iv,1,3},{jv,1,3}]//Simplify 
MnuAZsqpolR=Sum[Rsqpol[[iv,jv]]//Simplify,{iv,1,3},{jv,1,3}]//Simplify 
MnuAZsqpol=MnuAZsqpolL+MnuAZsqpolR//Simplify

dsignudz= coefww Mnusqpol/.Union[ksub,gsub] //Simplify
dsignuAdz= coefww MnuAsqpol/.Union[ksub,gsub] //Simplify
dsignuAZdz= coefww MnuAZsqpol/.Union[ksub,gsub] //Simplify

dsdznupb=hbcsq dsignudz/.numsub;
dsdznuApb=hbcsq dsignuAdz/.numsub;
dsdznuAZpb=hbcsq dsignuAZdz/.numsub;
Integrate[dsdznupb/. s-> 190^2,{(Subscript[c, θ]),-1,1}]
Integrate[dsdznuApb/. s-> 190^2,{(Subscript[c, θ]),-1,1}]
Integrate[dsdznuAZpb/. s-> 190^2,{(Subscript[c, θ]),-1,1}]
Plot[{dsdznupb/. s-> 190^2,dsdznuApb/. s-> 190^2,dsdznuAZpb/. s-> 190^2},{(Subscript[c, θ]),-1,1},PlotStyle->{{Blue,Dashed,Thick},{Red,Dotted,Thick},{Black,Thick}}]

dsignuAZdzL= coefww MnuAZsqpolL/.Union[ksub,gsub] //Simplify;
dsignuAZdzR= coefww MnuAZsqpolR/.Union[ksub,gsub] //Simplify;
dsdznuAZpbL=hbcsq dsignuAZdzL/.numsub;
dsdznuAZpbR=hbcsq dsignuAZdzR/.numsub;
Integrate[dsdznuAZpbR/. s-> 190^2,{(Subscript[c, θ]),-1,1}]
Integrate[dsdznuAZpbL/. s-> 190^2,{(Subscript[c, θ]),-1,1}]
Plot[{dsdznuAZpbR/. s-> 190^2,dsdznuAZpbL/. s-> 190^2},{(Subscript[c, θ]),-1,1},PlotStyle->{{Blue,Dashed,Thick},{Red,Dotted,Thick}}]

Plot[{Integrate[dsdznupb/. s->srt^2,{(Subscript[c, θ]),-1,1}],Integrate[dsdznuApb/. s->srt^2,{(Subscript[c, θ]),-1,1}],Integrate[dsdznuAZpb/. s->srt^2,{(Subscript[c, θ]),-1,1}]},{srt,161,200},PlotStyle->{{Blue,Dashed,Thick},{Red,Dotted,Thick},{Black,Thick}}]

Good qualitative agreement with figure 7.30 (the latter includes radiative corrections)

t  decay, Problem 7.10

FCClearScalarProducts;
{Coefv=((Subscript[p, fv])/(16 π (Subscript[m, t]^2))),ScalarProduct[p1n,p1n]=(Subscript[m, t]^2),ScalarProduct[p2n,p2n]=0,ScalarProduct[p3n,p3n]=(Subscript[M, W]^2),ScalarProduct[p1n,p2n]=(Subscript[p, fv]) (Subscript[m, t]),ScalarProduct[p1n,p3n]=(Subscript[m, t]) en3,ScalarProduct[p2n,p3n]=(1/2) ((Subscript[m, t]^2)-(Subscript[M, W]^2)),ScalarProduct[stn,stn]=-1,ScalarProduct[stn,p1n]=0,ScalarProduct[stn,p2n]=-(Subscript[p, fv]) (Subscript[c, θ]),ScalarProduct[stn,p3n]=+(Subscript[p, fv]) (Subscript[c, θ]),(Subscript[p, fv])=(((Subscript[m, t]^2)-(Subscript[M, W]^2))/(2 (Subscript[m, t]))),en3=(((Subscript[m, t]^2)+(Subscript[M, W]^2))/(2 (Subscript[m, t])))};
{ep3n=Momentum[Polarization[p3n]],
ep3sn= Momentum[Polarization[p3n]//ComplexConjugate]};
Pair[ep3n,Momentum[p1n]]=Pair[ep3n,Momentum[p2n]]=A
Pair[ep3sn,Momentum[p1n]]=Pair[ep3sn,Momentum[p2n]]=Astar

Polarized t, sum over W spins.

Pts=(1+ GA[5].GS[stn])/2
M1=Spinor[p2n,0].wfbfL[μ].Pts.Spinor[p1n, (Subscript[m, t])]
M1s=ComplexConjugate[M1]/.μ->ν
(FermionSpinSum[PolarizationSum[μ,ν,p3n]M1 M1s //EpsEvaluate //Explicit]//DiracTrick)
Mbsq=%/.DiracTrace->TR  //Contract  //Simplify
dsdz=Coefv Mbsq
(*Agrees with 7.10*)

W polarization for unpolarized t

(*M1=- I g Spinor[p2n,0]. GS[ep3sn].GA[7].Spinor[p1n, (Subscript[m, t])] /Sqrt[2]*)
M1=- I g Spinor[p2n,0]. DiracSlash[Polarization[p3n]//ComplexConjugate].GA[7].Spinor[p1n, (Subscript[m, t])] /Sqrt[2]
M1s=ComplexConjugate[M1]
Mbsq=FermionSpinSum[M1 M1s  //Explicit]/2/.DiracTrace->TR   //Simplify
Mbsqr=Mbsq/.{LeviCivita[][p1n,p2n,ep3sn,ep3n]->crterm}

Mbsq0=Mbsqr /.{crterm->0,Astar->(Subscript[p, fv]) (Subscript[m, t])/(Subscript[M, W]),A->(Subscript[p, fv]) (Subscript[m, t])/(Subscript[M, W])}//Simplify
MbsqL=Mbsqr /.{crterm->-I (Subscript[m, t]) (Subscript[p, fv]),Astar->0,A->0}//Simplify
MbsqR=Mbsqr /.{crterm->I (Subscript[m, t])(Subscript[p, fv]),Astar->0,A->0}//Simplify
F0=Mbsq0/(Mbsq0+MbsqL+MbsqR)//Simplify

Z pole asymmetries for massless fermions
Use general kinematic setup for 2 by 2 scattering

FCClearScalarProducts;
SetMandelstam[s, t, u, pg1,pg2, -pg3, -pg4,mg1,mg2,mg3,mg4];
{kgf=kfc[s,mg3,mg4],kgi=kfc[s,mg1,mg2],eg1=e1fc[s,mg1,mg2],eg2=e1fc[s,mg2,mg1],eg3=e1fc[s,mg3,mg4],eg4=e1fc[s,mg4,mg3],βg1=kgi/eg1,βg2=kgi/eg2, βg3=kgf/eg3, βg4=kgf/eg4};
coefg=kgf/(32 Pi s kgi); (*coefficient for ((dσ)/dz) *)
cmvar={t->mg1^2+mg3^2-2 eg1 eg3 + 2 kgi kgf cθ, u->mg1^2+mg4^2-2 eg1 eg4 - 2 kgi kgf cθ};
{thresholdi={s->(mg1+mg2)^2},thresholdf={s->(mg3+mg4)^2}};
ctheta = {cθ->Cos[(Subscript[θ, CM])]};

{mg1=0,mg2=0,mg3=0,mg4=0};  (*define masses*)

FB asymmetry

ampfb=Spinor[-pg2].zfbfg[μ,gV,gA].Spinor[pg1] (Spinor[pg3].zfbfg[μ,gV,gA].Spinor[-pg4] /.intofin);
sigfb=FermionSpinSum[ampfb(ComplexConjugate[ampfb]/.{μ->ν})]/.DiracTrace->TR// Contract //Simplify;
sigmaF=Integrate[sigfb/.cmvar, {cθ,0,1}];
%/.display
sigmaB=Integrate[sigfb/.cmvar, {cθ,-1,0}];
%/.display
AFB=(sigmaF-sigmaB)/(sigmaF+sigmaB)//Simplify  (*(7.172)*);
%/.display
AFB/.asympar;
%/.display

initial polarization

amph[h_]:=Spinor[-pg2].zfbfg[μ,gV,gA].GA[h].Spinor[pg1] (Spinor[pg3].zfbfg[μ,gV,gA].Spinor[-pg4] /.intofin)
sighf[h_]:=FermionSpinSum[amph[h] (ComplexConjugate[amph[h]]/.{μ->ν})]/.DiracTrace->TR// Contract //Simplify;

sigmaL=Integrate[sighf[7]/.cmvar, {cθ,-1,1}];
%/.display
sigmaR=Integrate[sighf[6]/.cmvar, {cθ,-1,1}];
%/.display
ALR=((sigmaL-sigmaR)/(sigmaL+sigmaR)//Simplify);  (*(7.175)*)
%/.display
ALR/.asympar;
%/.display

FBLR

sigmaLF=Integrate[sighf[7]/.cmvar, {cθ,0,1}];
%/.display
sigmaLB=Integrate[sighf[7]/.cmvar, {cθ,-1,0}];
%/.display
sigmaRF=Integrate[sighf[6]/.cmvar, {cθ,0,1}];
%/.display
sigmaRB=Integrate[sighf[6]/.cmvar, {cθ,-1,0}];
%/.display
AFBLR=(sigmaLF-sigmaLB-sigmaRF+sigmaRB)/(sigmaLF+sigmaLB+sigmaRF+sigmaRB) //Simplify (*(7.176)*);
%/.display
AFBLR/.asympar;
%/.display

final polarization

ampf[tauhel_]:=Spinor[-pg2].zfbfg[μ,gV,gA].Spinor[pg1] (Spinor[pg3].zfbfg[μ,gV,gA].GA[tauhel].Spinor[-pg4] /.intofin)
sigf[tauhel_]:=FermionSpinSum[ampf[tauhel] (ComplexConjugate[ampf[tauhel]]/.{μ->ν})]/.DiracTrace->TR// Contract //Simplify;

sigmafL=Integrate[sigf[7]/.cmvar, {cθ,-1,1}];
%/.display
sigmafR=Integrate[sigf[6]/.cmvar, {cθ,-1,1}];
%/.display
Ptau=((sigmafR-sigmafL)/(sigmafR+sigmafL)//Simplify);
%/.display    (*averaged over angle (7.173)*)
Ptau/.asympar;
%/.display

sigtR=((sigf[6]/.cmvar)/.cθ->z)//Simplify;
sigtL=((sigf[7]/.cmvar)/.cθ->z)//Simplify;
Ptauz=((sigtR-sigtL)/(sigtR+sigtL)//Simplify);
%/.display

Ptauzr=(Ptauz/.(z^2)+1->z/R)//FullSimplify;
%/.display
Ptauzr//.asympar//FullSimplify;
%/.display

TOC