NLO corrections to the R-ratio

This notebook contains the evaluation of the real and virtual corrections to the R-ration. The virtual diagrams are evaluated using package X, see https://packagex.hepforge.org

T. Becher, Nov. 2020

Auxiliary routines to perform the ϵ-expansion of the phase-space integrals

Distributions

DistApply[ex_,u_]:=Collect[ex,{DiracDelta[u],PlusDist[f_,u]}]//.distRules[u]

distRules[z_]:={a_ DiracDelta[z]:>(a/.z->0),a_. PlusDist[f_,z]:>f(a-(a/.z->0))};

nmax=3;

distExp[z_]:=(z^(-1+a_))->(1/a)DiracDelta[z]+Sum[((a^n)/(n!))PlusDist[(((Log[z])^n)/z),z],{n,0,nmax}];

Example

exampl=(z^(-1+ϵ))(((1-z))^ϵ);

Here we can just do the integration

fullRes=Integrate[exampl,{z,0,1},Assumptions->Re[ϵ]>0]

Expand in the integrand using distribution identity:

exampl2=exampl/.distExp[z]

exampl3=Series[exampl2,{ϵ,0,2}]//Normal

exampl4=DistApply[exampl3,z]

example5=Integrate[exampl4,{z,0,1}]

Series[fullRes,{ϵ,0,2}]//FullSimplify//Normal

Real emission corrections

Parameterization of the 3-particle phase space

This has to be integrated over the three-momenta (Subscript[p, 1]) and (Subscript[p, 2])

phase=((((1/(((2π))^(d-1)))))^3)(1/(8 E1 E2 E3))DiracDelta[Q-E1-E2-E3](((2π))^d);

solidAngle[d_]=((2(π^(d/2)))/(Gamma[d/2]));

Only one nontrivial angular integral: the angle between the three-vectors  (Subscript[p, 1]) and (Subscript[p, 2])

phase2=solidAng[d-1] solidAng[d-2] (E1^(d-2))(E2^(d-2))((Sin[θ])^(d-3)) phase/.E3->(Sqrt[(E1^2)+(E2^2)+2 E1 E2 Cos[θ]])

sol=Solve[0==E1+E2-Q+(Sqrt[(E1^2)+(E2^2)+2 E1 E2 Cos[θ]])/.Cos[θ]->Cθ,Cθ][[1]]//Simplify

fprime=Simplify[D[E1+E2-Q+(Sqrt[(E1^2)+(E2^2)+2 E1 E2 Cos[θ]])/.Cos[θ]->Cθ,Cθ]/.sol,{Q>E1+E2}]

Let us rewrite the integral as an integral over the Cθ = Cos[θ]; note  dθ ((Sin(θ))^a)== dcos(θ)  ((Sin(θ))^(a-1))=dcos(θ)  ((((Sin^2)(θ)))^((a-1)/2))

phase3=FullSimplify[(1/fprime)phase2/. ((Sin[θ])^(a_))->(((1-((Cθ)^2)))^((a-1)/2))/.Cos[θ]->Cθ/.sol,{Q>E1+E2,E1>0,E2>0}]

phase4=(2^(4-3 d)) (π^(3-2 d)) ((((2 E1+2 E2-Q) Q (-2 E1+Q) (-2 E2+Q)))^((1/2) (-4+d))) solidAng[-2+d] solidAng[-1+d]//Simplify

FullSimplify[phase4/phase3//Together//Factor,{Q>E1+E2,E1>0,E2>0}]

varch=Solve[{y1==1-((2E1)/Q),y2==1-((2E2)/Q)}][[1]]//Factor

jac=Det[{D[{E1,E2}/.varch,y1],D[{E1,E2}/.varch,y2]}]

phase5=FullSimplify[jac phase4/.varch/.d->4-2ϵ,{Q>0,y1>0,y2>0,y1<1,y2<1}]

phase6=FullSimplify[phase5/. -1+y1+y2->-y3,{Q>0}]//PowerExpand

phase7=phase6/.solidAng[a_]->solidAngle[a]//FullSimplify

prefps3=((Q^2)/(128(π^3)))(((((4π)/(Q^2))))^(2ϵ)) ;

phase8=phase7/prefps3//PowerExpand//Simplify

Amplitude and phase space

3-particle phase space, parameterized in terms of (Subscript[y, i]) = 1- ((2(Subscript[E, i]))/Q). 
Note that energy conservation gives (Subscript[y, 1])+(Subscript[y, 2])+(Subscript[y, 3]) = 1.

The integral goes over the region 0 ≤(Subscript[y, 1])+(Subscript[y, 2]) ≤1.

PS3=(Q2/(128 (π^3)))(((((4π)/Q2)))^(2ϵ))(1/(Gamma[2-2ϵ]))(((y1 y2 y3))^(-ϵ));

We want to normalize to the tree-level in d dimensions. To do so we also need the two-body phase space (Exercise)

PS2=(((2^(-4+4 ϵ)) (π^(-(1/2)+ϵ)) (Q2^(-ϵ)))/(Gamma[(3/2)-ϵ]));

qqbargAmpSquared=qqbarAmpSquared ((CF (g^2)((μ2)^ϵ))/Q2) ((4 (y3-y1 y2 ϵ)+2(1-ϵ)((y1^2)+(y2^2)))/(y1 y2));

prefact=((αs CF)/(4π))(((((4π μ2)/Q2)))^ϵ)(1/(Gamma[1-ϵ]));

To get the cross section in units of the LO cross section, we compute:

rest=(1/prefact)((PS3/PS2)//FullSimplify//PowerExpand) (qqbargAmpSquared/qqbarAmpSquared)/.g->(Sqrt[4π αs])//PowerExpand

... this needs to be integrated over (Subscript[y, 1])+(Subscript[y, 2]) ≤ 1.

Integrate[(y1^(-1-ϵ)),{y1,0,1}]

Phase-space integration: distribution expansion

rest2=rest/.y3->1-y2-y1//PowerExpand

Variable change: (Subscript[y, 2])->(1-(Subscript[y, 1])) u. Both variables u and y1 run from 0 ... 1

rest3=(1-y1) rest2/.y2->(1-y1) u/.(((u (1-y1)))^(-1-ϵ))->(u^(-1-ϵ)) (((1-y1))^(-1-ϵ))

Singularities at u = 0 and (Subscript[y, 1]) = 0, but the two integrations are decoupled.

rest4=Series[rest3/.distExp[u]/.distExp[y1],{ϵ,0,0}]//Normal

rest5=DistApply[DistApply[rest4,u],y1]

rest6=Integrate[rest5,{y1,0,1}]

rest7=Integrate[rest6,{u,0,1}]

rest8=Series[((E^(EulerGamma ϵ))/(Gamma[1-ϵ])) rest7,{ϵ,0,0}]//Normal//Expand

The final result is

σqqbg=σqqbar*((CF αs )/(4 π))(((((4 π (E^(-EulerGamma)) (µ^2))/Q2)))^ϵ)(19+(4/(ϵ^2))+(6/ϵ)-((7 (π^2))/3));

Alternative: subtraction scheme

If (Subscript[y, i]) -> 0 then the (Subscript[E, i])= Q/2. This is only possible if the other two particles are flying collinearly in the other direction.

If  (Subscript[y, 1]) -> 0 AND  (Subscript[y, 2]) -> 0 then (Subscript[E, 3]) -> 0, the gluon becomes soft.

rest

collinear13=Series[rest/.y3->1-y1-y2,{y2,0,0}]//Normal

collinear23=Series[rest/.y3->1-y1-y2,{y1,0,0}]//Normal

soft3=Series[collinear13,{y1,0,0}]//Normal

soft3alt=Series[collinear23,{y2,0,0}]//Normal

These are the regions in which the cross section becomes singular. Subtracting them, we get a finite cross section which can be integrated numerically:

finite=Series[rest-collinear13-collinear23+soft3/.y3->1-y1-y2,{ϵ,0,0}]

.... at higher orders in ϵ there are contributions.

Note: need to add the soft part because is contained twice in the collinear results.

Let us integrate the divergent pieces

collInt13=Integrate[Integrate[collinear13,{y2,0,1-y1},Assumptions->Re[ϵ]<0],{y1,0,1},Assumptions->Re[ϵ]<0]

softInt=Integrate[Integrate[soft3,{y2,0,1-y1},Assumptions->Re[ϵ]<0],{y1,0,1},Assumptions->Re[ϵ]<0]

sum=Series[2*collInt13-softInt,{ϵ,0,0}]//Normal

This reproduces the integral from above!

Loop corrections to (SuperPlus[e])(SuperMinus[e])⟶ q (OverBar[q])

Needs["X`"]

Package X

External line corrections: fermion self energy

Numerator of the diagram

self=CF (g^(2 ))DiracMatrix[(Subscript[γ, μ]),γ.(p+k),(Subscript[γ, ν])](-(Subscript[𝕘, μ,ν])k.k +ξ (Subscript[k, μ]) (Subscript[k, ν]));

(1/((((k.k))^2)(((p+k))^2)))

self2=self//Contract

self3=self2//FermionLineExpand

Add denominators and loop integration

self4=LoopIntegrate[self3,k,{k,0,2},{k+p,0}]

self5=LoopRefine[self4]

self4onShell=self4/.p.p->0

self5OnShell=LoopRefine[self4onShell]

Note that the package suppresses a factor

extra=((I (E^(-EulerGamma ϵ)))/(((4π))^(d/2)));

Vertex correction

triangle=(I^6)CF (g^(2 )) ⟨𝓊[p2,0],(Subscript[γ, α]),γ.(-p2+k),(Subscript[γ, μ]),γ.(p1+k),(Subscript[γ, β]),𝓊[p1,0]⟩(-(Subscript[𝕘, α,β])k.k +ξ (Subscript[k, α]) (Subscript[k, β]))

triangle2=(1/(CF (g^2)))triangle//Contract

triangle3=LoopIntegrate[triangle2,k,{k,0,2},{k+p1,0},{k-p2,0}]/.p1.p1->0/.p2.p2->0

triangle4=LoopRefine[triangle3]

triangle5=(1/(⟨𝓊[p2,0],(Subscript[γ, μ]),𝓊[p1,0]⟩ ))triangle4/. p1.p2->Q2/2

Note: Log[-((µ^2)/Q2)] has an imaginary part. The  iϵ ‘s in the propagators lead to Q2 ≡ Q2 + iϵ :

Therefore   Log[-((µ^2)/Q2)] = Log[((µ^2)/Q2)]+ i π

triangle6=triangle5/.Log[-((µ^2)/Q2)]->Log[((µ^2)/Q2)]+I π//Expand//Simplify

triangle7=(((((µ^2)/Q2)))^ϵ)(Series[(((((µ^2)/Q2)))^(-ϵ))triangle6,{ϵ,0,0}]//Simplify//Normal)

extra=((I (E^(-EulerGamma ϵ)))/(((4π))^(d/2)));

And the Feynman rules give i * Amp and we divided by (g^2) CF

prefact=(1/I)extra (g^2)CF/.g->(Sqrt[4π αs])/.d->4-2ϵ//Simplify

So the full result is

triangle8=((CF αs)/(4π))(((((4π (E^(-EulerGamma))(µ^2))/Q2)))^ϵ) (-8-(2/(ϵ^2))+((-3-2 I π)/ϵ)-3 I π+((7 (π^2))/6)) ⟨𝓊[p2,0],(Subscript[γ, μ]),𝓊[p1,0]⟩;

Including the electrons, squaring and writing the result in terms of the tree-level amplitude in d dimensions, we have:

σvirtual = ((CF αs)/(4π))HoldForm[(((((4π (E^(-EulerGamma))(µ^2))/Q2)))^ϵ) ]σ0 (2*(-8-(2/(ϵ^2))+((-3)/ϵ)+((7 (π^2))/6))//Expand)

(19+(4/(ϵ^2))+(6/ϵ)-((7 (π^2))/3))+(-16-(4/(ϵ^2))-(6/ϵ)+((7 (π^2))/3))//Expand

Checks

Test

integ=(u^(-1-ϵ)) (((1-y1))^(-ϵ)) (((1-u (1-y1)-y1))^(-ϵ)) (y1^(-1-ϵ)) (4 (1-u (1-y1)-y1-ϵ u (1-y1) y1)+2 (1-ϵ) ((u^2) (((1-y1))^2)+(y1^2)));

Simpler

integ=(u^(-1-ϵ))  (y1^(-1-ϵ)) (4 (1-u (1-y1)-y1-ϵ u (1-y1) y1)+2 (1-ϵ) ((u^2) (((1-y1))^2)+(y1^2)))

integ=(u^(-1-ϵ))  (y1^(-1-ϵ)) (4 (1-u (1-y1)-y1));

integ=(u^(-1-ϵ))  (y1^(-1-ϵ)) (1-u -y1);

NIntegrate[integ/.ϵ->-0.467,{y1,0,1},{u,0,1},PrecisionGoal->3]

numer=(integ/((u^(-1-ϵ)) (y1^(-1-ϵ)) ))

subtru=numer-(numer/.u->0)

sub=(numer/.u->0)Integrate[(u^(-1-ϵ)),{u,0,1},Assumptions->Re[ϵ]<0]

Full result is

NIntegrate[(u^(-1-ϵ)) (y1^(-1-ϵ))subtru/.ϵ->-0.467,{y1,0,1},{u,0,1},PrecisionGoal->3]+NIntegrate[(y1^(-1-ϵ))sub/.ϵ->-0.467,{y1,0,1},PrecisionGoal->3]

subtry1=subtru-(subtru/.y1->0)

intSubSub=Integrate[Series[(u^(-1-ϵ)) (y1^(-1-ϵ))subtry1,{ϵ,0,0}]//Normal,{u,0,1},{y1,0,1}]

Full result is

(subtru/.y1->0)

Integrate[(y1^(-1-ϵ)),{y1,0,1},Assumptions->Re[ϵ]<0]

NIntegrate[(u^(-1-ϵ)) (y1^(-1-ϵ))subtry1/.ϵ->-0.467,{y1,0,1},{u,0,1},PrecisionGoal->3]+(Integrate[(y1^(-1-ϵ)),{y1,0,1},Assumptions->Re[ϵ]<0]/.ϵ->-0.467) NIntegrate[(u^(-1-ϵ)) (subtru/.y1->0)/.ϵ->-0.467,{y1,0,1},{u,0,1},PrecisionGoal->3]+NIntegrate[(y1^(-1-ϵ))sub/.ϵ->-0.467,{y1,0,1},PrecisionGoal->3]

subtrY=(Integrate[(y1^(-1-ϵ)),{y1,0,1},Assumptions->Re[ϵ]<0]) Integrate[Series[(u^(-1-ϵ)) (subtru/.y1->0),{ϵ,0,1}]//Normal,{u,0,1}]

subInt=Series[Integrate[(y1^(-1-ϵ)) sub,{y1,0,1},Assumptions->Re[ϵ]<0],{ϵ,0,0}]//Normal

Add everything up...

subInt+subtrY+intSubSub//Expand

Test

integ

integ2=Series[(u^(-1-ϵ)) (1-u-y1) (y1^(-1-ϵ))/.distExp[u]/.distExp[y1],{ϵ,0,0}]

distExp

test=Integrate[integ,{y1,0,1},Assumptions->Re[ϵ]<0]

test2=Integrate[test,{u,0,1},Assumptions->Re[ϵ]<0]

Series[test2,{ϵ,0,0}]

integ2=Series[(1/(y1^(-1-ϵ)))integ/.distExp[u],{ϵ,0,1}]//Normal

one=(y1^(-1-ϵ))(DistCollect[integ2,u]/.distRules[u])/.distExp[y1]

two=DistCollect[one,y1]/.distRules[y1];

three=Integrate[two,{u,0,1}]

four=Integrate[three,{y1,0,1}]//Expand

Checks

ampAlt=(((x1^2)+(x2^2)-ϵ (x3^2))/((1-x1)(1-x2)))/.x1->1-y1/.x2->1-y2/.x3->1-y3/.y3->1-y1-y2//Simplify

((4 (y3-y1 y2 ϵ)+2(1-ϵ)((y1^2)+(y2^2)))/(y1 y2))-2ampAlt/.y3->1-y1-y2//Simplify

qqbargAmpSquaredPaper=((2CF (gs^2) ((μ2)^ϵ))/Q2)(((x1^2)+(x2^2)-ϵ (x3^2))/((1-x1)(1-x2)))/.x1->1-y1/.x2->1-y2/.x3->1-y3/.y3->1-y1-y2//Simplify

Übung

aha=Integrate[(x^(-1+ϵ))(((1-x))^(2+ϵ)),{x,0,1},GenerateConditions->False]

bebe=Integrate[((Sin[ϕ])^(-1+ϵ)),{ϕ,0,π},GenerateConditions->False]

Series[bebe,{ϵ,0,0}]//FullSimplify

bebe=Integrate[((Sin[a x])^(-1+ϵ)),{x,0,1},GenerateConditions->False]