Local Unitarity for (SuperPlus[e]) (SuperMinus[e]) -> (SuperStar[γ])-> d (OverBar[d]) at NLO

Global definitions

Resources

Let the integrator always display error message as it provides an estimate of the integral error.

Off[General::stop];

Dot products

LDot[k_,l_]:=k[[1]]l[[1]]-k[[2]]l[[2]]-k[[3]]l[[3]]-k[[4]]l[[4]]

VDot[k_,l_]:=k[[1]]l[[1]]+k[[2]]l[[2]]+k[[3]]l[[3]]

MyPrint[msg_]:=Print[Style[msg,FontFamily->"Monaco"]]

Energy[k_List]:=k[[1]]
Xcomp[k_List]:=k[[2]]
Ycomp[k_List]:=k[[3]]
Zcomp[k_List]:=k[[4]]

GeVtoPbConversion=0.389379304 (10^9);

SP[x_List,y_List]:=LDot[x,y]

Convolution function h(t)

h[t_,σ_]:=Exp[-((t^2)/(σ^2))](2/(Sqrt[π]σ))

Test it :

Integrate[h[t,12],{t,0,∞}]

Parametrisations

SphericalMap[x1_,x2_,x3_]:=Module[
{r=(x1/(1-x1)),θ=x2 π,ϕ=x3 2π,y,drOverdx=(1/(((x1-1))^2))},
<|"vector"->{
r Sin[θ]Cos[ϕ],
r Sin[θ]Sin[ϕ],
r Cos[θ]
},
"jacobian"->(π) (2π)drOverdx (r^2)Sin[θ]
|>
]

CartesianMap[x1_,x2_,x3_]:=Module[
{Map=(y/(1-y))-((1-y)/y),dypOvery=(1/(y^2))+(1/(((y-1))^2))},
<|"vector"->{
Map/.{y->x1},
Map/.{y->x2},
Map/.{y->x3}
},
"jacobian"->(dypOvery/.{y->x1})(dypOvery/.{y->x2})(dypOvery/.{y->x3})
|>
]

Test it:

(*Simplify[NIntegrate[If[
(((SphericalMap[x1,x2,x3][["vector"]][[1]])^2)
+((SphericalMap[x1,x2,x3][["vector"]][[2]])^2)
+((SphericalMap[x1,x2,x3][["vector"]][[3]])^2))<4,SphericalMap[x1,x2,x3][["jacobian"]],0],
{x1,0,1},{x2,0,1},{x3,0,1}]-(4/3)π 8]*)

(*Simplify[NIntegrate[If[
(((CartesianMap[x1,x2,x3][["vector"]][[1]])^2)
+((CartesianMap[x1,x2,x3][["vector"]][[2]])^2)
+((CartesianMap[x1,x2,x3][["vector"]][[3]])^2))<4,CartesianMap[x1,x2,x3][["jacobian"]],0],
{x1,0,1},{x2,0,1},{x3,0,1}]-(4/3)π 8]*)

Facility to manipulate tensor polynomials

(*RemoveZeroCoefficients[coefs_,threshold_]:=Module[{result},
result=DeleteCases[coefs,{___,sym[_]}/;True];
DeleteCases[result,x_/;Abs[x[[-1]][[1]]]≤threshold]
]*)

(*RemoveZeroCoefficients[coefs_]:=RemoveZeroCoefficients[coefs,0]*)

Facility to evaluate symmetric tensor polynomials

(*EvaluatePolynomial[coefs_,loopmomenta_]:=Module[{},
Total[Table[
Product[
Product[
loopmomenta[[IcoefsForOneLoop]][[onepower+1]]
,{onepower,term[[IcoefsForOneLoop]]/.{sym[x_]:>x}}]
,{IcoefsForOneLoop,Length[term[[;;-2]]]}]*term[[-1]][[1]]
,{term,coefs}
]]
]*)

Facility to convert a polynomial into a list of symmetric coefficients

Parameters

For clarity we regroup here all couplings:

AllNLOCouplings=(gs^2)(ge^4)(Qd^2)(Qe^2);

SelectedP1={1,0,0,1};

SelectedP2={1,0,0,-1};

NumericalValues={
ge->N[Sqrt[4π (1/(132.507`16))]],
gs->N[Sqrt[4π 0.118`16]],
Qe->1,Qd->(1/3),d->4,Nc->3,TF->(1/2),
p1->SelectedP1,
p2->SelectedP2
};

Leading Order

We are interested in the one supergraph contributing to (SuperPlus[e]) (SuperMinus[e]) -> (SuperStar[γ])-> d (OverBar[d]) at leading order. 
We will work in the rest frame where of the off-shell photon, that is (Q^μ)=((Q^0),0,0,0), (Subscript[p, 1]^μ)=(((Q^0)/2),0,0,((Q^0)/2)) and (Subscript[p, 1]^μ)=(((Q^0)/2),0,0,-((Q^0)/2)). 
For convenience, we will set (Q^0)=2 from now on

Polynomial Numerator

LONumeratorPolynomial=(k0^2*(64*ge^4*Nc*Qd^2*Qe^2-32*d*ge^4*Nc*Qd^2*Qe^2)+kx^2*(-128*ge^4*Nc*Qd^2*Qe^2+32*d*ge^4*Nc*Qd^2*Qe^2)+ky^2*(-128*ge^4*Nc*Qd^2*Qe^2+32*d*ge^4*Nc*Qd^2*Qe^2)+kz^2*(-64*ge^4*Nc*Qd^2*Qe^2+32*d*ge^4*Nc*Qd^2*Qe^2)+k0*(-128*ge^4*Nc*Qd^2*Qe^2+64*d*ge^4*Nc*Qd^2*Qe^2))/.NumericalValues

LOIntegrand[k_List]:=(GeVtoPbConversion((((1/(-I))))^2)(1/(((2π))^4))(1/(2SP[p1+p2,p1+p2]))(((-2π I))/(2Abs[Energy[k]]))(((-2π I))/(2Abs[(Energy[p1]+Energy[p2])-Energy[k]]))
((1/((2 2)))((((1/(SP[p1+p2,p1+p2]))))^2))LONumeratorPolynomial/.{k0->k[[1]],kx->k[[2]],ky->k[[3]],kz->k[[4]]})/.NumericalValues

LOIntegrand[{1,2,3,4}]

Defining the complete integrand

Now the procedure for solving the rescaling delta

LORescalingSolutionParametric=Solve[
FullSimplify[Sqrt[(t^2)INkdotk]+Sqrt[(t^2)INkdotk]-INQ0==0,
Assumptions->{t>0,INkdotk>0,INQ0>0}],{t}
];
LORescalingSolutionParametricJacobian=
Abs[(D[FullSimplify[Sqrt[(t^2)INkdotk]+Sqrt[(t^2)INkdotk]-INQ0,
Assumptions->{t>0,INkdotk>0,INQ0>0}],t])/.{t->LORescalingSolutionParametric}];
LORescalingSolution[Q0_,kvec_]:=(LORescalingSolutionParametric/.{INkdotk->VDot[kvec,kvec],INQ0->Q0})
LORescalingSolutionJacobian[Q0_,kvec_]:=Abs[(LORescalingSolutionParametricJacobian/.{INkdotk->VDot[kvec,kvec],INQ0->Q0})]

We can now proceed with the integration

LOIntegrandFunction[xs_,OptionsPattern[DEBUG->False]]:=Module[
{
kmap=SphericalMap[xs[[1]],xs[[2]],xs[[3]]],
jacobian=1,Integrand=1,bubble,kvec,k0,Q0,kmu,QMinuskmu,l0,rescaling,rescalingJacobian,
δ=0,λ=1,σ,DEBUG=OptionValue[DEBUG],p1=SelectedP1,p2=SelectedP2
},
kvec=kmap["vector"];
jacobian=jacobian*kmap["jacobian"];
Q0=Energy[p1]+Energy[p2];
σ=Q0;
rescaling=(t/.LORescalingSolution[Q0,kvec][[1]]);
If[DEBUG,MyPrint["rescaling="<>ToString[rescaling,InputForm]];];
rescalingJacobian=LORescalingSolutionJacobian[Q0,kvec];
kvec=rescaling*kvec;
k0=Sqrt[VDot[kvec,kvec]];
kmu={k0,kvec[[1]],kvec[[2]],kvec[[3]]};
QMinuskmu={Q0-k0,-kvec[[1]],-kvec[[2]],-kvec[[3]]};
If[DEBUG,MyPrint["kmu=("<>ToString[k0,InputForm]<>","<>ToString[kvec,InputForm]<>")"];];
If[DEBUG,MyPrint["|k|+|k|="<>ToString[Sqrt[VDot[kvec,kvec]]+Sqrt[VDot[kvec,kvec]],InputForm]<>" vs Q0="<>ToString[Q0,InputForm]]];
If[DEBUG,MyPrint["|k|="<>ToString[Sqrt[VDot[kvec,kvec]]]<>" vs k0="<>ToString[k0,InputForm]]];
If[DEBUG,MyPrint["|k|="<>ToString[Sqrt[VDot[kvec,kvec]]]<>" vs Q0-k0="<>ToString[Q0-k0,InputForm]]];
Integrand=Integrand*LOIntegrand[kmu];
If[DEBUG,MyPrint["LO Integrand="<>ToString[LOIntegrand[kmu],InputForm]];];
Integrand=Integrand*h[rescaling,σ]*((rescaling^3));
If[DEBUG,MyPrint["h[rescaling,σ]="<>ToString[h[rescaling,σ],InputForm]];];
If[DEBUG,MyPrint["(rescaling^3)="<>ToString[(rescaling^3),InputForm]];];
Integrand=Integrand*(1/rescalingJacobian);
If[DEBUG,MyPrint["1/(rescaling jacobian)="<>ToString[1/rescalingJacobian,InputForm]];];
Integrand=Integrand*jacobian;
If[DEBUG,MyPrint["jacobian="<>ToString[jacobian,InputForm]];];
If[DEBUG,MyPrint["Integrand="<>ToString[Integrand,InputForm]];];
Re[Integrand]
]

Test the integrand now

LOIntegrandFunction[{0.4,0.2,0.3},DEBUG->True]

Performing the numerical integration

ClearAll[ThisIntegrand,CutFunction];
ThisIntegrand[x1,x2,x3]=LOIntegrandFunction[{x1,x2,x3},DEBUG->False];
(*MCPoints={};*)
rBoundMin=0.0+(10^(-5));rBoundMax=1.0-(10^(-5));
BoundMin=0.0+(10^(-5));BoundMax=1.0-(10^(-5));
IntermediateResult=0;
LOResults=Association[Table[
IntermediateResult=NIntegrate[
ThisIntegrand[x1,x2,x3],
{x1,rBoundMin,rBoundMax},{x2,BoundMin,BoundMax},{x3,BoundMin,BoundMax},
Method->"AdaptiveMonteCarlo",PrecisionGoal->4,MaxPoints->statistics,MaxRecursion->1000
(*,EvaluationMonitor:>AppendTo[MCPoints,{{x1,x2,x3},RealCutRescalingIntegrand[{x1,x2,x3},0,2,False],VirtualCutRescalingIntegrand[{x1,x2,x3},0,2,False]}]*)
];
MyPrint["Result for "<>ToString[statistics]<>" sample points: "<>ToString[IntermediateResult]];
statistics->IntermediateResult
,{statistics,{1000,5000,10000,50000,100000,500000,1000000}}]
]

FinalLOResult = LOResults[Max[Keys[LOResults]]]

Cross-section :  ( 0.7744 +- 0.0004323 pb )*100

LOResultBenchmark=7744.0;

LOResultRelDiff=Abs[(FinalLOResult-LOResultBenchmark)/LOResultBenchmark]*100.0

Next-to-Leading Order

We are interested in the two supergraphs contributing to (SuperPlus[e]) (SuperMinus[e]) -> (SuperStar[γ])-> d (OverBar[d]) at next-to-leading order.
We will work in the frame (Subscript[p, 1]^μ)=(1,0,0,1) and (Subscript[p, 1]^μ)=(1,0,0,-1).

Nested self-energy

Real-emission self-energy (RBubble)

Polynomial Numerator

Note below that we used below that the polarisation sum of gluons if -(g^(μν)) which would not reproduce the actual squared diagram in MadGraph for example, but it is ok thanks to the ward identity

RBubbleNumeratorPolynomial=((8*I)*(-2+d)*ge^4*gs^2*((-2+d)*k0^3*l0+4*ky^2*l0-2*d*ky^2*l0+4*kz^2*l0-2*d*kz^2*l0+(-4+d)*kx^3*lx+kx*((-4+d)*ky^2+d*kz^2)*lx-4*ky^3*ly+d*ky^3*ly+d*ky*kz^2*ly-k0*((-6+d)*kx^2*l0+(-6+d)*ky^2*l0-4*(-2+d)*kx*lx-4*(-2+d)*ky*ly+(-2+d)*kz*(kz*l0-4*lz))-6*ky^2*kz*lz+d*ky^2*kz*lz-2*kz^3*lz+d*kz^3*lz+kx^2*(-2*(-2+d)*l0+(-4+d)*ky*ly+(-6+d)*kz*lz)-k0^2*(2*(-2+d)*l0-2*kz*lz+d*(kx*lx+ky*ly+kz*lz)))*(-1+Nc^2)*Qd^2*Qe^2*TF)/.NumericalValues;

RBubbleME[k_List,l_List]:=RBubbleNumeratorPolynomial/.{
k0->k[[1]],kx->k[[2]],ky->k[[3]],kz->k[[4]],
l0->l[[1]],lx->l[[2]],ly->l[[3]],lz->l[[4]]}

RBubbleIntegrand[k_List,l_List]:=Module[{},
(GeVtoPbConversion(1/(((2π))^8))(1/(2SP[p1+p2,p1+p2]))(((-2π I))/(2Abs[Energy[l]]))(((-2π I))/(2Abs[Energy[k]-Energy[l]]))(((-2π I))/(2Abs[(Energy[p1]+Energy[p2])-Energy[k]]))(
((((1/(SP[p1+p2,p1+p2]))(1/(SP[k,k]))))^2)RBubbleME[k,l]
))/.NumericalValues
]

RBubbleIntegrand[{0.1,0.2,0.3,0.4},{0.4,0.5,0.6,0.7}]

Defining the complete integrand

RBubbleRescalingSolutionParametric=Solve[
FullSimplify[Sqrt[(t^2)INldotl]+Sqrt[(t^2)INlkdotlk]+Sqrt[(t^2)INkdotk]-INQ0==0,
Assumptions->{t>0,INkdotk>0,INlkdotlk>0,INldotl>0,INQ0>0}],{t}
];
RBubbleRescalingSolutionParametricJacobian=
Abs[(D[FullSimplify[Sqrt[(t^2)INldotl]+Sqrt[(t^2)INlkdotlk]+Sqrt[(t^2)INkdotk]-INQ0,
Assumptions->{t>0,INkdotk>0,INlkdotlk>0,INldotl>0,INQ0>0}],t])/.{t->RBubbleRescalingSolutionParametric}];
RBubbleRescalingSolution[Q0_,kvec_,lvec_]:=(RBubbleRescalingSolutionParametric/.{INldotl->VDot[lvec,lvec],INlkdotlk->VDot[kvec-lvec,kvec-lvec],INkdotk->VDot[-kvec,-kvec],INQ0->Q0})
RBubbleRescalingSolutionJacobian[Q0_,kvec_,lvec_]:=(RBubbleRescalingSolutionParametricJacobian/.{INldotl->VDot[lvec,lvec],INlkdotlk->VDot[kvec-lvec,kvec-lvec],INkdotk->VDot[-kvec,-kvec],INQ0->Q0})

RBubbleIntegrandFunction[xs_,OptionsPattern[{DEBUG->False,SHIFTED->False}]]:=Module[
{
kmap=SphericalMap[xs[[1]],xs[[2]],xs[[3]]],
(*kmap=CartesianMap[xs[[1]],xs[[2]],xs[[3]]],*)
lmap=SphericalMap[xs[[4]],xs[[5]],xs[[6]]],
(*lmap=CartesianMap[xs[[4]],xs[[5]],xs[[6]]],*)
jacobian=1,Integrand=1,kvec,lvec,kmu,lmu,k0,l0,QMinuskmu,kmlmu,
rescaling,rescalingJacobian,Q0,σ,
DEBUG=OptionValue[DEBUG],SHIFTED=OptionValue[SHIFTED],
p1=SelectedP1,p2=SelectedP2
},
Q0=Energy[p1]+Energy[p2];
σ=Q0;
If[DEBUG,MyPrint["xs="<>ToString[xs]];];
lvec=lmap["vector"];
jacobian=jacobian*lmap["jacobian"];
kvec=kmap["vector"];
(* Apply shift if asked for *)
lvec=If[SHIFTED,lvec+(1/2)kvec,lvec];
jacobian=jacobian*kmap["jacobian"];
rescaling=(t/.RBubbleRescalingSolution[Q0,lvec,kvec][[1]]);
If[DEBUG,MyPrint["rescaling="<>ToString[rescaling]];];
rescalingJacobian=RBubbleRescalingSolutionJacobian[Q0,lvec,kvec];
lvec=rescaling*lvec;
kvec=rescaling*kvec;
l0=Sqrt[VDot[lvec,lvec]];
(*Must pick negative solution*)
k0=-Sqrt[VDot[kvec,kvec]]+Q0;
kmu={k0,kvec[[1]],kvec[[2]],kvec[[3]]};
QMinuskmu={Q0-k0,-kvec[[1]],-kvec[[2]],-kvec[[3]]};
lmu={l0,lvec[[1]],lvec[[2]],lvec[[3]]};
kmlmu=kmu-lmu;
If[DEBUG,MyPrint["lmu=("<>ToString[l0]<>","<>ToString[lvec]<>")"];];
If[DEBUG,MyPrint["kmu=("<>ToString[k0]<>","<>ToString[kvec]<>")"];];
If[DEBUG,MyPrint["-lmu+kmu=("<>ToString[-l0+k0]<>","<>ToString[-lvec+kvec]<>")"];];
If[DEBUG,MyPrint["|l|+|l-k|+|k|="<>ToString[Sqrt[VDot[lvec,lvec]]+Sqrt[VDot[lvec-kvec,lvec-kvec]]+Sqrt[VDot[kvec,kvec]]]<>" vs Q0="<>ToString[Q0]]];
If[DEBUG,MyPrint["|l-k|="<>ToString[Sqrt[VDot[lvec-kvec,lvec-kvec]]]<>" vs -(l0-k0)="<>ToString[-(l0-k0)]]];
If[DEBUG,MyPrint["|l|="<>ToString[Sqrt[VDot[lvec,lvec]]]<>" vs l0="<>ToString[l0]]];
If[DEBUG,MyPrint["|k|="<>ToString[Sqrt[VDot[kvec,kvec]]]<>" vs k0="<>ToString[k0]]];
If[DEBUG,MyPrint["|k|="<>ToString[Sqrt[VDot[kvec,kvec]]]<>" vs -(k0-Q0)="<>ToString[-(k0-Q0)]]];
Integrand=Integrand*RBubbleIntegrand[kmu,lmu];
If[DEBUG,MyPrint["(1/Eb1)*(1/Eb2)="<>ToString[(1/(Abs[Energy[lmu]]))(1/(Abs[Energy[kmu]-Energy[lmu]])),InputForm]];];
If[DEBUG,MyPrint["Integrand="<>ToString[Integrand]];];
If[DEBUG,MyPrint["Integrand*jac="<>ToString[Integrand*jacobian,InputForm]];];
Integrand=Integrand*h[rescaling,σ]*((rescaling^6));
If[DEBUG,MyPrint["h[rescaling,σ]="<>ToString[h[rescaling,σ]]];];
If[DEBUG,MyPrint["(rescaling^6)="<>ToString[(rescaling^6)]];];
Integrand=Integrand*(1/rescalingJacobian);
If[DEBUG,MyPrint["1/(rescaling jacobian)="<>ToString[1/rescalingJacobian]];];
Integrand=Integrand*jacobian;
If[DEBUG,MyPrint["jacobian="<>ToString[jacobian]];];
If[DEBUG,MyPrint["Integrand="<>ToString[Integrand]];];
Re[Integrand]
]

RBubbleIntegrandFunction[{0.1,0.2,0.3,0.4,0.5,0.6},DEBUG->True]

Virtual on-shell self-energy (VBubble)

Polynomial Numerator

VBubblePolynomialSelfEnergyToTheLeft=(((-4*I)*(-2+d)*ge^4*gs^2*(k0+l0-m0)*((-2+d)*(-2+k0)*l0^3-4*ky*lx^2*ly+d*ky*lx^2*ly-4*ky*ly^3+d*ky*ly^3-2*kz*lx^2*lz+d*kz*lx^2*lz-2*kz*ly^2*lz+d*kz*ly^2*lz-4*ky*ly*lz^2+d*ky*ly*lz^2-2*kz*lz^3+d*kz*lz^3+8*lx^2*m0-4*d*lx^2*m0-4*k0*lx^2*m0+2*d*k0*lx^2*m0+8*ly^2*m0-4*d*ly^2*m0-4*k0*ly^2*m0+2*d*k0*ly^2*m0+8*lz^2*m0-4*d*lz^2*m0-4*k0*lz^2*m0+2*d*k0*lz^2*m0-4*ky*ly*m0^2+d*ky*ly*m0^2-2*kz*lz*m0^2+d*kz*lz*m0^2-8*lx*m0*mx+4*d*lx*m0*mx+4*k0*lx*m0*mx-2*d*k0*lx*m0*mx+4*ky*ly*mx^2-d*ky*ly*mx^2+2*kz*lz*mx^2-d*kz*lz*mx^2+8*ky*lx^2*my-2*d*ky*lx^2*my+8*ky*ly^2*my-2*d*ky*ly^2*my+8*ky*lz^2*my-2*d*ky*lz^2*my-8*ly*m0*my+4*d*ly*m0*my+4*k0*ly*m0*my-2*d*k0*ly*m0*my-8*ky*lx*mx*my+2*d*ky*lx*mx*my-4*ky*ly*my^2+d*ky*ly*my^2+2*kz*lz*my^2-d*kz*lz*my^2-l0^2*((-4+d)*kx*(lx-2*mx)+(-4+d)*ky*(ly-2*my)+(-2+d)*(2*(-2+k0)*m0+kz*(lz-2*mz)))+4*kz*lx^2*mz-2*d*kz*lx^2*mz+4*kz*ly^2*mz-2*d*kz*ly^2*mz+4*kz*lz^2*mz-2*d*kz*lz^2*mz-8*lz*m0*mz+4*d*lz*m0*mz+4*k0*lz*m0*mz-2*d*k0*lz*m0*mz-4*kz*lx*mx*mz+2*d*kz*lx*mx*mz-4*kz*ly*my*mz+2*d*kz*ly*my*mz-8*ky*lz*my*mz+2*d*ky*lz*my*mz+4*ky*ly*mz^2-d*ky*ly*mz^2-2*kz*lz*mz^2+d*kz*lz*mz^2+l0*(-((-2+d)*(-2+k0)*lx^2)-(-2+d)*(-2+k0)*ly^2-4*lz^2+2*d*lz^2+2*k0*lz^2-d*k0*lz^2+4*m0^2-2*d*m0^2-2*k0*m0^2+d*k0*m0^2+8*kx*m0*mx-2*d*kx*m0*mx+4*mx^2-2*d*mx^2-2*k0*mx^2+d*k0*mx^2+8*ky*m0*my-2*d*ky*m0*my+4*my^2-2*d*my^2-2*k0*my^2+d*k0*my^2+4*kz*m0*mz-2*d*kz*m0*mz+4*mz^2-2*d*mz^2-2*k0*mz^2+d*k0*mz^2)+(-4+d)*kx*(lx^3-2*lx^2*mx-2*mx*(ly^2-ly*my+lz*(lz-mz))+lx*(ly^2+lz^2+m0^2+mx^2-my^2-mz^2)))*(-1+Nc^2)*Qd^2*Qe^2*TF)/Pi)/.NumericalValues

VBubblePolynomialSelfEnergyToTheRight=(((4*I)*(-2+d)*ge^4*gs^2*(k0+l0-m0)*((-2+d)*(-2+k0)*l0^3-4*ky*lx^2*ly+d*ky*lx^2*ly-4*ky*ly^3+d*ky*ly^3-2*kz*lx^2*lz+d*kz*lx^2*lz-2*kz*ly^2*lz+d*kz*ly^2*lz-4*ky*ly*lz^2+d*ky*ly*lz^2-2*kz*lz^3+d*kz*lz^3+8*lx^2*m0-4*d*lx^2*m0-4*k0*lx^2*m0+2*d*k0*lx^2*m0+8*ly^2*m0-4*d*ly^2*m0-4*k0*ly^2*m0+2*d*k0*ly^2*m0+8*lz^2*m0-4*d*lz^2*m0-4*k0*lz^2*m0+2*d*k0*lz^2*m0-4*ky*ly*m0^2+d*ky*ly*m0^2-2*kz*lz*m0^2+d*kz*lz*m0^2-8*lx*m0*mx+4*d*lx*m0*mx+4*k0*lx*m0*mx-2*d*k0*lx*m0*mx+4*ky*ly*mx^2-d*ky*ly*mx^2+2*kz*lz*mx^2-d*kz*lz*mx^2+8*ky*lx^2*my-2*d*ky*lx^2*my+8*ky*ly^2*my-2*d*ky*ly^2*my+8*ky*lz^2*my-2*d*ky*lz^2*my-8*ly*m0*my+4*d*ly*m0*my+4*k0*ly*m0*my-2*d*k0*ly*m0*my-8*ky*lx*mx*my+2*d*ky*lx*mx*my-4*ky*ly*my^2+d*ky*ly*my^2+2*kz*lz*my^2-d*kz*lz*my^2-l0^2*((-4+d)*kx*(lx-2*mx)+(-4+d)*ky*(ly-2*my)+(-2+d)*(2*(-2+k0)*m0+kz*(lz-2*mz)))+4*kz*lx^2*mz-2*d*kz*lx^2*mz+4*kz*ly^2*mz-2*d*kz*ly^2*mz+4*kz*lz^2*mz-2*d*kz*lz^2*mz-8*lz*m0*mz+4*d*lz*m0*mz+4*k0*lz*m0*mz-2*d*k0*lz*m0*mz-4*kz*lx*mx*mz+2*d*kz*lx*mx*mz-4*kz*ly*my*mz+2*d*kz*ly*my*mz-8*ky*lz*my*mz+2*d*ky*lz*my*mz+4*ky*ly*mz^2-d*ky*ly*mz^2-2*kz*lz*mz^2+d*kz*lz*mz^2+l0*(-((-2+d)*(-2+k0)*lx^2)-(-2+d)*(-2+k0)*ly^2-4*lz^2+2*d*lz^2+2*k0*lz^2-d*k0*lz^2+4*m0^2-2*d*m0^2-2*k0*m0^2+d*k0*m0^2+8*kx*m0*mx-2*d*kx*m0*mx+4*mx^2-2*d*mx^2-2*k0*mx^2+d*k0*mx^2+8*ky*m0*my-2*d*ky*m0*my+4*my^2-2*d*my^2-2*k0*my^2+d*k0*my^2+4*kz*m0*mz-2*d*kz*m0*mz+4*mz^2-2*d*mz^2-2*k0*mz^2+d*k0*mz^2)+(-4+d)*kx*(lx^3-2*lx^2*mx-2*mx*(ly^2-ly*my+lz*(lz-mz))+lx*(ly^2+lz^2+m0^2+mx^2-my^2-mz^2)))*(-1+Nc^2)*Qd^2*Qe^2*TF)/Pi)/.NumericalValues

Note: the UV counterterm considered integrates directly into the correct finite part to reproduce the MSbar and on-shell renormalisation scheme

VBubblePolynomialSelfEnergyUVCT=(-256*ge^4*gs^2*(k0^4-2*k0^3*(1+2*l0)+k0^2*(-2*kz^2+4*(-8+2*l0+l0^2)+2*kz*lz)+8*k0*(8+kz^2*l0-l0^2-kz*l0*lz)+8*kz*(-(kz*(-4+l0^2))+l0^2*lz))*Pi^2*Qd^2*Qe^2)/.NumericalValues

Now we can build the complete matrix element contributions from the above numerators, and test them with some predefined input

kvecTest={0.1,0.2,0.3};
kTest={Sqrt[VDot[kvecTest,kvecTest]],kvecTest[[1]],kvecTest[[2]],kvecTest[[3]]};
lvecTest={0.4,0.5,0.6};
lTest={Sqrt[VDot[lvecTest,lvecTest]],lvecTest[[1]],lvecTest[[2]],lvecTest[[3]]};

VBubbleLeftIntegrand[kIN_,lIN_]:=Block[
{
k,l,p1,p2,qE,lvec,kvec,kE,zstar,lE,mE,m,mvec
},
k=kIN;
l=lIN;
p1=SelectedP1;
p2=SelectedP2;
qE=Energy[k];
lvec={Xcomp[l],Ycomp[l],Zcomp[l]};
kvec={Xcomp[k],Ycomp[k],Zcomp[k]};
lE=Sqrt[VDot[lvec,lvec]];
kE=Sqrt[VDot[lvec-kvec,lvec-kvec]]+lE;
mvec=lvec-kvec;
mE=-Sqrt[VDot[mvec,mvec]];
m={mE,mvec[[1]],mvec[[2]],mvec[[3]]};
zstar=kE;
GeVtoPbConversion((((1/(-I))))^1)(1/(((2π))^4))(1/(2LDot[p1+p2,p1+p2]))(((-2π I))/(2Abs[(Energy[p1]+Energy[p2])-Energy[k]]))(
((1/(((2π))^4))(((-2π I))/(2qE))((((-2π I))/(2lE))(((-2π I))/(2(kE-lE))))(1/(((zstar^2)-(qE^2))))((1/(SP[p1+p2,p1+p2])))(1/(((zstar^2)-(qE^2))))((1/(SP[p1+p2,p1+p2]))))
(
VBubblePolynomialSelfEnergyToTheLeft/.{
k0->qE,kx->kvec[[1]],ky->kvec[[2]],kz->kvec[[3]],
l0->lE,lx->lvec[[1]],ly->lvec[[2]],lz->lvec[[3]],
m0->mE,mx->mvec[[1]],my->mvec[[2]],mz->mvec[[3]]
}
)
)
]

VBubbleLeftIntegrand[kTest,lTest]

VBubbleRightIntegrand[kIN_,lIN_]:=Block[
{
k,l,p1,p2,qE,lvec,kvec,kE,zstar,lE,mE,m,mvec
},
k=kIN;
l=lIN;
p1=SelectedP1;
p2=SelectedP2;
qE=Energy[k];
lvec={Xcomp[l],Ycomp[l],Zcomp[l]};
kvec={Xcomp[k],Ycomp[k],Zcomp[k]};
lE=Sqrt[VDot[lvec,lvec]];
kE=Sqrt[VDot[lvec-kvec,lvec-kvec]]+lE;
mvec=lvec-kvec;
mE=-Sqrt[VDot[mvec,mvec]];
m={mE,mvec[[1]],mvec[[2]],mvec[[3]]};
zstar=kE;
GeVtoPbConversion((((1/(-I))))^1)(1/(((2π))^4))(1/(2LDot[p1+p2,p1+p2]))(((-2π I))/(2Abs[(Energy[p1]+Energy[p2])-Energy[k]]))(
((1/(((2π))^4))(((-2π I))/(2qE))((((-2π I))/(2lE))(((-2π I))/(2(kE-lE))))(1/(((zstar^2)-(qE^2))))((1/(SP[p1+p2,p1+p2])))(1/(((zstar^2)-(qE^2))))((1/(SP[p1+p2,p1+p2]))))
(
VBubblePolynomialSelfEnergyToTheRight/.{
k0->qE,kx->kvec[[1]],ky->kvec[[2]],kz->kvec[[3]],
l0->-lE,lx->lvec[[1]],ly->lvec[[2]],lz->lvec[[3]],
m0->-mE,mx->mvec[[1]],my->mvec[[2]],mz->mvec[[3]]
}
)
)
]

VBubbleRightIntegrand[kTest,lTest]

VBubbleUVCTEnergyIntegrated[k_List,l_List]:=Module[
{
lmu,lvec,kvec,kmu,
Omega,lSopervec,MUV,p1,p2
},
p1=SelectedP1;
p2=SelectedP2;
MUV=2 (Energy[p1]+Energy[p2]);
lmu={l0,Xcomp[l],Ycomp[l],Zcomp[l]};
lvec={Xcomp[l],Ycomp[l],Zcomp[l]};
kmu={Energy[k],Xcomp[k],Ycomp[k],Zcomp[k]};
kvec={Xcomp[k],Ycomp[k],Zcomp[k]};
lSopervec=lvec-(1/2)kvec;
Omega=Sqrt[VDot[lSopervec,lSopervec]+(MUV^2)];
(* Notice that the factor 2π arising from the integration over the loop momenta energy and which needed to be added "manually" in LxB4DUVEnergyIntegrated[...] is in this case already part of the definition of BubbleUVCTCoefsValues.
Ben, you'll have to decide if you want it in or out, but below I don't multiply by 2π since it is already in BubbleUVCTCoefsValues. *)
((1/((3-1)!))D[(VBubblePolynomialSelfEnergyUVCT/.{
l0->(lSoper0+(1/2)Energy[k]),lx->lvec[[1]], ly->lvec[[2]], lz->lvec[[3]],
k0->Energy[k],kx->kvec[[1]], ky->kvec[[2]], kz->kvec[[3]]
})(1/(((lSoper0+Omega))^3)),{lSoper0,3-1}])/.{lSoper0->Omega}
]

VBubbleUVCTIntegrand[kIN_List,lIN_List]:=Block[
{
k,l,p1,p2,q0,lvec,kvec,k0,l0
},
k=kIN;
l=lIN;
p1=SelectedP1;
p2=SelectedP2;
qE=Energy[k];
GeVtoPbConversion((((1/(-I))))^1)(1/(((2π))^4))(1/(2LDot[p1+p2,p1+p2]))(((-2π I))/(2Abs[(Energy[p1]+Energy[p2])-Energy[k]]))(
(((1/(SP[p1+p2,p1+p2])))(1/(((2π))^4))(((-2π I))/(2qE))(1/(2π I))((1/(SP[p1+p2,p1+p2]))))
(
VBubbleUVCTEnergyIntegrated[k,l]
)
)
]

VBubbleUVCTIntegrand[kTest,lTest]

Now we can combine together all pieces

alaLTDSoperLRDifferentiatedVBubbleIntegrand[k_,l_]:={VBubbleLeftIntegrand[k,l],VBubbleRightIntegrand[k,l],VBubbleUVCTIntegrand[k,l]}

Defining the complete integrand

VBubbleRescalingSolutionParametric=Solve[
FullSimplify[Sqrt[(t^2)INkdotk]+Sqrt[(t^2)INkdotk]-INQ0==0,
Assumptions->{t>0,INkdotk>0,INQ0>0}],{t}
];
VBubbleRescalingSolutionParametricJacobian=
Abs[(D[FullSimplify[Sqrt[(t^2)INkdotk]+Sqrt[(t^2)INkdotk]-INQ0,
Assumptions->{t>0,INkdotk>0,INQ0>0}],t])/.{t->VBubbleRescalingSolutionParametric}];
VBubbleRescalingSolution[Q0_,kvec_]:=(VBubbleRescalingSolutionParametric/.{INkdotk->VDot[kvec,kvec],INQ0->Q0})
VBubbleRescalingSolutionJacobian[Q0_,kvec_]:=Abs[(VBubbleRescalingSolutionParametricJacobian/.{INkdotk->VDot[kvec,kvec],INQ0->Q0})]

alaLTDSoperLRDifferentiatedVBubbleIntegrandFunction[xs_,OptionsPattern[{DEBUG->False,SHIFTED->False,SHIFTEDUV->False}]]:=Module[
{
kmap=SphericalMap[xs[[1]],xs[[2]],xs[[3]]],
(*kmap=CartesianMap[xs[[1]],xs[[2]],xs[[3]]],*)
lmap=SphericalMap[xs[[4]],xs[[5]],xs[[6]]],
(*lmap=CartesianMap[xs[[4]],xs[[5]],xs[[6]]],*)
jacobian=1,Integrand=1,kvec,lvec,kmu,lmu,k0,l0,rescaling,rescalingJacobian,Q0,σ,VBubbleIntegrandEval,
p1=SelectedP1,p2=SelectedP2,
DEBUG=OptionValue[DEBUG],SHIFTED=OptionValue[SHIFTED],SHIFTEDUV=OptionValue[SHIFTEDUV]
},
Q0=Energy[p1]+Energy[p2];
σ=Q0;
If[DEBUG,MyPrint["xs="<>ToString[xs]];];
lvec=lmap["vector"];
jacobian=jacobian*lmap["jacobian"];
kvec=kmap["vector"];
lvec=If[SHIFTED,lvec+(1/2)kvec,lvec];
jacobian=jacobian*kmap["jacobian"];
rescaling=(t/.VBubbleRescalingSolution[Q0,kvec][[1]]);
If[DEBUG,MyPrint["rescaling="<>ToString[rescaling]];];
rescalingJacobian=VBubbleRescalingSolutionJacobian[Q0,kvec];
lvec=rescaling*lvec;
kvec=rescaling*kvec;
(* Set l0 to infinity to make sure it is not used. *)
l0=∞;
(*Must pick negative solution*)
k0=-Sqrt[VDot[kvec,kvec]]+Q0;
kmu={k0,kvec[[1]],kvec[[2]],kvec[[3]]};
lmu={l0,lvec[[1]],lvec[[2]],lvec[[3]]};
If[DEBUG,MyPrint["lmu=("<>ToString[l0]<>","<>ToString[lvec]<>")"];];
If[DEBUG,MyPrint["kmu=("<>ToString[k0]<>","<>ToString[kvec]<>")"];];
If[DEBUG,MyPrint["-lmu+kmu=("<>ToString[-l0+k0]<>","<>ToString[-lvec+kvec]<>")"];];
If[DEBUG,MyPrint["|k|+|k|="<>ToString[Sqrt[VDot[kvec,kvec]]+Sqrt[VDot[kvec,kvec]]]<>" vs Q0="<>ToString[Q0]]];
If[DEBUG,MyPrint["|k|="<>ToString[Sqrt[VDot[kvec,kvec]]]<>" vs k0="<>ToString[k0]]];
If[DEBUG,MyPrint["|k|="<>ToString[Sqrt[VDot[kvec,kvec]]]<>" vs Q0-k0="<>ToString[Q0-k0]]];
VBubbleIntegrandEval=alaLTDSoperLRDifferentiatedVBubbleIntegrand[kmu,lmu];
Integrand=Integrand*Total[VBubbleIntegrandEval];
If[DEBUG,MyPrint["Integrand alaLTDVBubbleIntegrand*jac  ="<>ToString[VBubbleIntegrandEval[[1]]*jacobian,InputForm]];];
If[DEBUG,MyPrint["UV CT alaLTDVBubbleIntegrand*jac          ="<>ToString[VBubbleIntegrandEval[[3]]*jacobian,InputForm]];];
If[And[DEBUG,VBubbleIntegrandEval[[3]]≠0],MyPrint["Ratio                        ="<>ToString[(VBubbleIntegrandEval[[1]]+VBubbleIntegrandEval[[2]])/VBubbleIntegrandEval[[3]],InputForm]];];
Integrand=Integrand*h[rescaling,σ]*((rescaling^6));
If[DEBUG,MyPrint["h[rescaling,σ]="<>ToString[h[rescaling,σ]]];];
If[DEBUG,MyPrint["(rescaling^6)="<>ToString[(rescaling^6)]];];
Integrand=Integrand*(1/rescalingJacobian);
If[DEBUG,MyPrint["1/(rescaling jacobian)="<>ToString[1/rescalingJacobian]];];
Integrand=Integrand*jacobian;
If[DEBUG,MyPrint["jacobian="<>ToString[jacobian]];];
If[DEBUG,MyPrint["Integrand="<>ToString[Integrand]];];
Re[Integrand]
]

alaLTDSoperLRDifferentiatedVBubbleIntegrandFunction[{0.6,0.4,0.3,0.33,0.40001,0.30001},DEBUG->True]

Numerical integration of the self-energy supergraph

ThisIntegrand[x1,x2,x3,x4,x5,x6]=RBubbleIntegrandFunction[{x1,x2,x3,x4,x5,x6},DEBUG->False]+alaLTDSoperLRDifferentiatedVBubbleIntegrandFunction[{x1,x2,x3,x4,x5,x6},DEBUG->False];
(*MCPoints={};*)
rBoundMin=0.0;rBoundMax=1.0;
BoundMin=0.0;BoundMax=1.0;
IntermediateResult=0;
FinalNLOVBubbleResult=Association[Table[
IntermediateResult=NIntegrate[
ThisIntegrand[x1,x2,x3,x4,x5,x6],
{x1,rBoundMin,rBoundMax},{x2,BoundMin,BoundMax},{x3,BoundMin,BoundMax},
{x4,rBoundMin,rBoundMax},{x5,BoundMin,BoundMax},{x6,BoundMin,BoundMax},
Method->"AdaptiveMonteCarlo",PrecisionGoal->5,AccuracyGoal->5,MaxPoints->statistics,MaxRecursion->10000
(*,EvaluationMonitor:>AppendTo[MCPoints,{{x1,x2,x3,x4,x5,x6},RealCutRescalingIntegrand[{x1,x2,x3,x4,x5,x6},0,2,False],VirtualCutRescalingIntegrand[{x1,x2,x3,x4,x5,x6},0,2,False]}]*)
];
Print["Result for "<>ToString[statistics]<>" sample points: "<>ToString[IntermediateResult]];
statistics->IntermediateResult
,{statistics, {10000,100000}}]
]

FinalNLOVBubbleResult= NLOVBubbleResults[Max[Keys[NLOVBubbleResults]]]

Cross-section : -360.732 +- 0.451046 [pb]

NLOVBubbleResultBenchmark=-360.0;

NLOVBubbleResultRelDiff=Abs[(FinalNLOVBubbleResult-NLOVBubbleResultBenchmark)/NLOVBubbleResultBenchmark]*100.0

Double Triangle

Virtual L x B* (LxB)

Polynomial Numerator

LxBNumeratorPolynomial=((-16*I)*ge^4*gs^2*(8*ky^2*l0-6*d*ky^2*l0+d^2*ky^2*l0+4*kz^2*l0-4*d*kz^2*l0+d^2*kz^2*l0+2*ky^2*l0^2-d*ky^2*l0^2+2*kz^2*l0^2-d*kz^2*l0^2-32*ky*ly+16*d*ky*ly-2*d^2*ky*ly-4*ky*l0*ly+2*d*ky*l0*ly-4*ky^2*ly^2+2*d*ky^2*ly^2-24*kz*lz+16*d*kz*lz-2*d^2*kz*lz-4*ky*kz*ly*lz+2*d*ky*kz*ly*lz-2*ky^2*lz^2+d*ky^2*lz^2-2*kz^2*lz^2+d*kz^2*lz^2+(-2+d)*k0^2*(-((-2+d)*l0)+l0^2-lz^2)+(-2+d)*kx^2*((-4+d)*l0-l0^2+2*lx^2+lz^2)+2*(-2+d)*k0*(-l0^2+kx*lx+ky*ly+l0*(-2+d-kx*lx-ky*ly)+lz^2)-2*kx*lx*(d^2+2*(8+l0+2*ky*ly+kz*lz)-d*(8+l0+2*ky*ly+kz*lz)))*(-1+Nc^2)*Qd^2*Qe^2*TF)/.NumericalValues

LxBME[k_,l_]:=Module[
{MENum,p1=SelectedP1,p2=SelectedP2 },
MENum =((((1/(SP[p1+p2,p1+p2]))))^2)LxBNumeratorPolynomial/.{
k0->k[[1]],kx->k[[2]],ky->k[[3]],kz->k[[4]],
l0->l[[1]],lx->l[[2]],ly->l[[3]],lz->l[[4]]
};
MENum
]

LxBMELTD[k_,l_]:=Module[
{LTD1,LTD2,LTD3,UV,Q0,MUV, kvec,lvec,k1,k2,k3,
p1=SelectedP1,p2=SelectedP2 
},
Q0=Energy[LTD1]+Energy[LTD2];
MUV=2 Q0;
lvec={Xcomp[l],Ycomp[l],Zcomp[l]};
kvec={Xcomp[k],Ycomp[k],Zcomp[k]};
k1={Sqrt[VDot[kvec,kvec]],Xcomp[k],Ycomp[k],Zcomp[k]};
LTD1=LxBME[k1,l](1/(2Sqrt[VDot[kvec,kvec]]))(1/(SP[l-k1,l-k1]))(1/(SP[k1-p1-p2,k1-p1-p2]));
k2={Sqrt[VDot[kvec,kvec]]+Energy[p1]+Energy[p2],Xcomp[k],Ycomp[k],Zcomp[k]};
LTD2=LxBME[k2,l](1/(SP[k2,k2]))(1/(SP[l-k2,l-k2]))(1/(2Sqrt[VDot[kvec,kvec]]));
k3={Sqrt[VDot[kvec-lvec,kvec-lvec]]+Energy[l],Xcomp[k],Ycomp[k],Zcomp[k]};
LTD3=LxBME[k3,l](1/(SP[k3,k3]))(1/(2Sqrt[VDot[kvec-lvec,kvec-lvec]]))(1/(SP[k3-p1-p2,k3-p1-p2]));
(1/(((2π))^4))(-2π I ){LTD1,LTD2,LTD3}
]

LxBIntegrand[k_,l_]:=Block[{p1=SelectedP1,p2=SelectedP2},
GeVtoPbConversion (1/(((2π))^4))(1/(2LDot[p1+p2,p1+p2]))(((-2π I))/(2Abs[Energy[l]]))(((-2π I))/(2Abs[(Energy[p1]+Energy[p2])-Energy[l]]))LxBMELTD[k,l]
]

LxBIntegrand[{0.1,0.2,0.3,0.4},{0.5,0.6,0.7,0.8}]

Now the 4D UV counterterm (whose integrated counterterpart is already such that it reproduces MSbar renormalisation conditions)

LxB4DUVNumeratorPolynomial=(32*ge^4*gs^2*((-32+k0^2-kx^2-ky^2-kz^2)*l0^2+2*kx^2*lx^2+4*kx*ky*lx*ly+2*ky^2*ly^2+2*k0*(kx*lx+ky*ly)+l0*(-64+2*k0^2-4*kx^2-4*ky^2-2*kz^2+d*(32-k0^2+kx^2+ky^2+kz^2)-2*k0*(kx*lx+ky*ly))+2*kx*kz*lx*lz+2*ky*kz*ly*lz+32*lz^2-k0^2*lz^2+kx^2*lz^2+ky^2*lz^2+kz^2*lz^2)*(-1+Nc^2)*Qd^2*Qe^2*TF)/.NumericalValues

LxB4DUVEnergyIntegrated[k_,l_]:=Module[
{
p1=SelectedP1,
p2=SelectedP2,
lmu={Energy[l],Xcomp[l],Ycomp[l],Zcomp[l]},
kvec={Xcomp[k],Ycomp[k],Zcomp[k]},
kmu={k0,Xcomp[k],Ycomp[k],Zcomp[k]
}
},(2π((1/((3-1)!))D[((LxB4DUVNumeratorPolynomial/.{kx->kvec[[1]],ky->kvec[[2]],kz->kvec[[3]],l0->lmu[[1]],lx->lmu[[2]],ly->lmu[[3]],lz->lmu[[4]]})/(((k0+Sqrt[VDot[kvec,kvec]+(MUV^2)]))^3)),{k0,3-1}])/.{k0->Sqrt[VDot[kvec,kvec]+(MUV^2)]})/.{MUV->2 (Energy[p1]+Energy[p2])}
]

LxB4DUVIntegrand[k_List,l_List]:=Module[
{p1=SelectedP1,p2=SelectedP2},
(GeVtoPbConversion (((1/(LDot[p1+p2,p1+p2])))
LxB4DUVEnergyIntegrated[k,l]
((1/(LDot[p1+p2,p1+p2])))(1/(((2π))^8))(1/(2LDot[p1+p2,p1+p2]))(((-2π I))/(2Abs[Energy[l]]))(((-2π I))/(2Abs[(Energy[p1]+Energy[p2])-Energy[l]]))))
]

LxB4DUVIntegrand[{0.1,0.2,0.3,0.4},{0.5,0.6,0.7,0.8}]

Now combine it with the UV

LxBCombinedIntegrand[k_,l_]:=Join[LxBIntegrand[k,l],{LxB4DUVIntegrand[k,l]}]

Defining the complete integrand

LxBRescalingSolutionParametric=Solve[
FullSimplify[Sqrt[(t^2)INldotl]+Sqrt[(t^2)INldotl]-INQ0==0,
Assumptions->{t>0,INldotl>0,INQ0>0}],{t}
];
LxBRescalingSolutionParametricJacobian=
Abs[(D[FullSimplify[Sqrt[(t^2)INldotl]+Sqrt[(t^2)INldotl]-INQ0,
Assumptions->{t>0,INldotl>0,INQ0>0}],t])/.{t->LxBRescalingSolutionParametric}];
LxBRescalingSolution[Q0_,lvec_]:=(LxBRescalingSolutionParametric/.{INldotl->VDot[lvec,lvec],INQ0->Q0})
LxBRescalingSolutionJacobian[Q0_,lvec_]:=Abs[(LxBRescalingSolutionParametricJacobian/.{INldotl->VDot[lvec,lvec],INQ0->Q0})]

LxBIntegrandFunction[xs_,OptionsPattern[{DEBUG->True}]]:=Module[
{
kmap=SphericalMap[xs[[1]],xs[[2]],xs[[3]]],
(*kmap=CartesianMap[xs[[1]],xs[[2]],xs[[3]]],*)
lmap=SphericalMap[xs[[4]],xs[[5]],xs[[6]]],
(*lmap=CartesianMap[xs[[4]],xs[[5]],xs[[6]]],*)
jacobian=1,Integrand=1,kvec,lvec,kmu,lmu,k0,l0,rescaling,rescalingJacobian,Q0,σ,LTDIntegrand,
p1=SelectedP1, p2=SelectedP2,
QMinuslmu,DEBUG=OptionValue[DEBUG]
},
Q0=Energy[p1]+Energy[p2];
σ=Q0;
If[DEBUG,MyPrint["xs="<>ToString[xs]];];
lvec=lmap["vector"];
jacobian=jacobian*lmap["jacobian"];
kvec=kmap["vector"];
jacobian=jacobian*kmap["jacobian"];
If[DEBUG,MyPrint["k jacobian="<>ToString[kmap["jacobian"],InputForm]];];
If[DEBUG,MyPrint["|k|="<>ToString[Sqrt[VDot[kvec,kvec]],InputForm]];];
rescaling=(t/.LxBRescalingSolution[Q0,lvec][[1]]);
If[DEBUG,MyPrint["rescaling="<>ToString[rescaling]];];
rescalingJacobian=LxBRescalingSolutionJacobian[Q0,lvec];
lvec=rescaling*lvec;
kvec=rescaling*kvec;
(* Set l0 to infinity to make sure it is not used. *)
k0=∞;
l0=Sqrt[VDot[lvec,lvec]];
kmu={k0,kvec[[1]],kvec[[2]],kvec[[3]]};
lmu={l0,lvec[[1]],lvec[[2]],lvec[[3]]};
QMinuslmu={Q0-l0,-lvec[[1]],-lvec[[2]],-lvec[[3]]};
If[DEBUG,MyPrint["lmu=("<>ToString[l0]<>","<>ToString[lvec]<>")"];];
If[DEBUG,MyPrint["kmu=("<>ToString[k0]<>","<>ToString[kvec]<>")"];];
If[DEBUG,MyPrint["-lmu+kmu=("<>ToString[-l0+k0]<>","<>ToString[-lvec+kvec]<>")"];];
If[DEBUG,MyPrint["|l|+|l|="<>ToString[Sqrt[VDot[lvec,lvec]]+Sqrt[VDot[lvec,lvec]]]<>" vs Q0="<>ToString[Q0]]];
If[DEBUG,MyPrint["|l|="<>ToString[Sqrt[VDot[lvec,lvec]]]<>" vs k0="<>ToString[l0]]];
If[DEBUG,MyPrint["|l|="<>ToString[Sqrt[VDot[lvec,lvec]]]<>" vs Q0-l0="<>ToString[Q0-l0]]];
LTDIntegrand=LxBCombinedIntegrand[kmu,lmu];
Integrand=Integrand*Total[LTDIntegrand];
If[DEBUG,MyPrint["Integrand LTD1  ="<>ToString[LTDIntegrand[[1]],InputForm]];];
If[DEBUG,MyPrint["Integrand LTD2  ="<>ToString[LTDIntegrand[[2]],InputForm]];];
If[DEBUG,MyPrint["Integrand LTD3  ="<>ToString[LTDIntegrand[[3]],InputForm]];];
If[DEBUG,MyPrint["Integrand UVCT  ="<>ToString[LTDIntegrand[[4]],InputForm]];];
If[DEBUG,MyPrint["Total[LTD]      ="<>ToString[Total[LTDIntegrand],InputForm]];];
If[DEBUG,MyPrint["Ratio LTD/UV    ="<>ToString[((Total[LTDIntegrand[[;;-2]]])/(LTDIntegrand[[4]])),InputForm]];];
If[DEBUG,MyPrint["Integrand LTD1*jac  ="<>ToString[LTDIntegrand[[1]]*jacobian,InputForm]];];
If[DEBUG,MyPrint["Integrand LTD2*jac  ="<>ToString[LTDIntegrand[[2]]*jacobian,InputForm]];];
If[DEBUG,MyPrint["Integrand LTD3*jac  ="<>ToString[LTDIntegrand[[3]]*jacobian,InputForm]];];
If[DEBUG,MyPrint["Integrand UVCT*jac  ="<>ToString[LTDIntegrand[[4]]*jacobian,InputForm]];];
If[DEBUG,MyPrint["Total[LTD]*jac      ="<>ToString[Total[LTDIntegrand]*jacobian,InputForm]];];
Integrand=Integrand*h[rescaling,σ]*((rescaling^6));
If[DEBUG,MyPrint["h[rescaling,σ]="<>ToString[h[rescaling,σ]]];];
If[DEBUG,MyPrint["(rescaling^6)="<>ToString[(rescaling^6)]];];
Integrand=Integrand*(1/rescalingJacobian);
If[DEBUG,MyPrint["1/(rescaling jacobian)="<>ToString[1/rescalingJacobian]];];
Integrand=Integrand*jacobian;
If[DEBUG,MyPrint["jacobian="<>ToString[jacobian]];];
If[DEBUG,MyPrint["Integrand="<>ToString[Integrand]];];
Re[Integrand]
]

LxBIntegrandFunction[{0.1,0.2,0.3,0.9999,0.5,0.6},DEBUG->True]

Virtual B x L* (BxL)

Polynomial Numerator

BxLNumeratorPolynomial=((16*I)*ge^4*gs^2*(2*(kx*lx+ky*ly)*(-d^2-2*(8+l0+kx*lx+ky*ly)+d*(8+l0+kx*lx+ky*ly))-2*(-2+d)*kz*(-6+d-kx*lx-ky*ly)*lz+(-2+d)*k0^2*(-2*l0+l0^2-lx^2-ly^2-lz^2)-(-2+d)*kz^2*(-2*l0+l0^2-lx^2-ly^2-lz^2)-(-2+d)*k0*((-2+d)*l0^2-2*kx*lx+4*lx^2-d*lx^2-2*ky*ly+4*ly^2-d*ly^2+l0*(4-2*d+2*kx*lx+2*ky*ly)+2*lz^2-d*lz^2))*(-1+Nc^2)*Qd^2*Qe^2*TF)/.NumericalValues

BxLME[k_,l_]:=Module[
{MENum,p1=SelectedP1,p2=SelectedP2 },
MENum =((((1/(SP[p1+p2,p1+p2]))))^2)BxLNumeratorPolynomial/.{
k0->k[[1]],kx->k[[2]],ky->k[[3]],kz->k[[4]],
l0->l[[1]],lx->l[[2]],ly->l[[3]],lz->l[[4]]
};
MENum
]

BxLMELTD[k_,l_]:=Module[
{LTD1,LTD2,LTD3,UV,kvec,lvec,l1,l2,l3,
p1=SelectedP1,p2=SelectedP2},
lvec={Xcomp[l],Ycomp[l],Zcomp[l]};
kvec={Xcomp[k],Ycomp[k],Zcomp[k]};
l1={Sqrt[VDot[lvec,lvec]],Xcomp[l],Ycomp[l],Zcomp[l]};
LTD1=BxLME[k,l1](1/(2Sqrt[VDot[lvec,lvec]]))(1/(SP[l1-k,l1-k]))(1/(SP[l1-p1-p2,l1-p1-p2]));
l2={Sqrt[VDot[lvec,lvec]]+Energy[p1]+Energy[p2],Xcomp[l],Ycomp[l],Zcomp[l]};
LTD2=BxLME[k,l2](1/(SP[l2,l2]))(1/(SP[l2-k,l2-k]))(1/(2Sqrt[VDot[lvec,lvec]]));
l3={Sqrt[VDot[lvec-kvec,lvec-kvec]]+Energy[k],Xcomp[l],Ycomp[l],Zcomp[l]};
LTD3=BxLME[k,l3](1/(SP[l3,l3]))(1/(2Sqrt[VDot[lvec-kvec,lvec-kvec]]))(1/(SP[l3-p1-p2,l3-p1-p2]));
(1/(((2π))^4))(+2π I ){LTD1,LTD2,LTD3}
]

BxLIntegrand[k_,l_]:=Block[{p1=SelectedP1,p2=SelectedP2},
GeVtoPbConversion(1/(((2π))^4))(1/(2LDot[p1+p2,p1+p2]))(((-2π I))/(2Abs[Energy[k]]))(((-2π I))/(2Abs[(Energy[p1]+Energy[p2])-Energy[k]]))BxLMELTD[k,l]
]

BxLIntegrand[{0.1,0.2,0.3,0.4},{0.5,0.6,0.7,0.8}]

Now the 4D UV counterterm (whose integrated counterterpart is already such that it reproduces MSbar renormalisation conditions)

BxL4DUVNumeratorPolynomial=(32*ge^4*gs^2*(2*(kx*lx+ky*ly)*(l0+kx*lx+ky*ly)+2*kz*(kx*lx+ky*ly)*lz+k0^2*(-32+l0^2-lx^2-ly^2-lz^2)+d*k0*(32-l0^2+lx^2+ly^2+lz^2)+kz^2*(32-l0^2+lx^2+ly^2+lz^2)-2*k0*(32-l0^2+kx*l0*lx+2*lx^2+ky*l0*ly+2*ly^2+lz^2))*(-1+Nc^2)*Qd^2*Qe^2*TF)/.NumericalValues

BxL4DUVEnergyIntegrated[k_,l_]:=Module[
{
p1=SelectedP1,
p2=SelectedP2,
kmu={Energy[k],Xcomp[k],Ycomp[k],Zcomp[k]},
lvec={Xcomp[l],Ycomp[l],Zcomp[l]},
lmu={l0,Xcomp[l],Ycomp[l],Zcomp[l]
}
},(2π((1/((3-1)!))D[((BxL4DUVNumeratorPolynomial/.{k0->kmu[[1]],kx->kmu[[2]],ky->kmu[[3]],kz->kmu[[4]],lx->lvec[[1]],ly->lvec[[2]],lz->lvec[[3]]})/(((l0+Sqrt[VDot[lvec,lvec]+(MUV^2)]))^3)),{l0,3-1}])/.{l0->Sqrt[VDot[lvec,lvec]+(MUV^2)]})/.{MUV->2 (Energy[p1]+Energy[p2])}
]

BxL4DUVIntegrand[k_List,l_List]:=Module[
{p1=SelectedP1,p2=SelectedP2},
(GeVtoPbConversion (((1/(LDot[p1+p2,p1+p2])))
BxL4DUVEnergyIntegrated[k,l]
((1/(LDot[p1+p2,p1+p2])))(1/(((2π))^8))(1/(2LDot[p1+p2,p1+p2]))(((-2π I))/(2Abs[Energy[k]]))(((-2π I))/(2Abs[(Energy[p1]+Energy[p2])-Energy[k]]))))
]

BxL4DUVIntegrand[{3,4,5,6},{7,8,9,10}]

Now combine it with the UV

BxLCombinedIntegrand[k_,l_]:=Join[BxLIntegrand[k,l],{BxL4DUVIntegrand[k,l]}]

Defining the complete integrand

BxLRescalingSolutionParametric=Solve[
FullSimplify[Sqrt[(t^2)INkdotk]+Sqrt[(t^2)INkdotk]-INQ0==0,
Assumptions->{t>0,INkdotk>0,INQ0>0}],{t}
];
BxLRescalingSolutionParametricJacobian=
Abs[(D[FullSimplify[Sqrt[(t^2)INkdotk]+Sqrt[(t^2)INkdotk]-INQ0,
Assumptions->{t>0,INkdotk>0,INQ0>0}],t])/.{t->BxLRescalingSolutionParametric}];
BxLRescalingSolution[Q0_,kvec_]:=(BxLRescalingSolutionParametric/.{INkdotk->VDot[kvec,kvec],INQ0->Q0})
BxLRescalingSolutionJacobian[Q0_,kvec_]:=Abs[(BxLRescalingSolutionParametricJacobian/.{INkdotk->VDot[kvec,kvec],INQ0->Q0})]

BxLIntegrandFunction[xs_,OptionsPattern[DEBUG->False]]:=Module[
{
kmap=SphericalMap[xs[[1]],xs[[2]],xs[[3]]],
(*kmap=CartesianMap[xs[[1]],xs[[2]],xs[[3]]],*)
lmap=SphericalMap[xs[[4]],xs[[5]],xs[[6]]],
(*lmap=CartesianMap[xs[[4]],xs[[5]],xs[[6]]],*)
jacobian=1,Integrand=1,kvec,lvec,kmu,lmu,k0,l0,rescaling,rescalingJacobian,Q0,σ,LTDIntegrand,
p1=SelectedP1,p2=SelectedP2,
DEBUG=OptionValue[DEBUG],QMinuskmu
},
Q0=Energy[p1]+Energy[p2];
σ=Q0;
If[DEBUG,MyPrint["xs="<>ToString[xs]];];
lvec=lmap["vector"];
jacobian=jacobian*lmap["jacobian"];
If[DEBUG,MyPrint["l jacobian="<>ToString[lmap["jacobian"],InputForm]];];
If[DEBUG,MyPrint["|l|="<>ToString[Sqrt[VDot[lvec,lvec]],InputForm]];];
kvec=kmap["vector"];
jacobian=jacobian*kmap["jacobian"];
rescaling=(t/.BxLRescalingSolution[Q0,kvec][[1]]);
If[DEBUG,MyPrint["rescaling="<>ToString[rescaling]];];
rescalingJacobian=BxLRescalingSolutionJacobian[Q0,kvec];
lvec=rescaling*lvec;
kvec=rescaling*kvec;
(* Set l0 to infinity to make sure it is not used. *)
l0=∞;
k0=Sqrt[VDot[kvec,kvec]];
kmu={k0,kvec[[1]],kvec[[2]],kvec[[3]]};
QMinuskmu={Q0-k0,-kvec[[1]],-kvec[[2]],-kvec[[3]]};
lmu={l0,lvec[[1]],lvec[[2]],lvec[[3]]};
If[DEBUG,MyPrint["lmu=("<>ToString[l0]<>","<>ToString[lvec]<>")"];];
If[DEBUG,MyPrint["kmu=("<>ToString[k0]<>","<>ToString[kvec]<>")"];];
If[DEBUG,MyPrint["-lmu+kmu=("<>ToString[-l0+k0]<>","<>ToString[-lvec+kvec]<>")"];];
If[DEBUG,MyPrint["|k|+|k|="<>ToString[Sqrt[VDot[kvec,kvec]]+Sqrt[VDot[kvec,kvec]]]<>" vs Q0="<>ToString[Q0]]];
If[DEBUG,MyPrint["|k|="<>ToString[Sqrt[VDot[kvec,kvec]]]<>" vs k0="<>ToString[k0]]];
If[DEBUG,MyPrint["|k|="<>ToString[Sqrt[VDot[kvec,kvec]]]<>" vs Q0-k0="<>ToString[Q0-k0]]];
LTDIntegrand=BxLCombinedIntegrand[kmu,lmu];
Integrand=Integrand*Total[LTDIntegrand];
If[DEBUG,MyPrint["Integrand LTD1  ="<>ToString[LTDIntegrand[[1]],InputForm]];];
If[DEBUG,MyPrint["Integrand LTD2  ="<>ToString[LTDIntegrand[[2]],InputForm]];];
If[DEBUG,MyPrint["Integrand LTD3  ="<>ToString[LTDIntegrand[[3]],InputForm]];];
If[DEBUG,MyPrint["Integrand UVCT  ="<>ToString[LTDIntegrand[[4]],InputForm]];];
If[DEBUG,MyPrint["Total[LTD]      ="<>ToString[Total[LTDIntegrand],InputForm]];];
If[DEBUG,MyPrint["Ratio LTD/UV    ="<>ToString[((Total[LTDIntegrand[[;;-2]]])/(LTDIntegrand[[4]])),InputForm]];];
If[DEBUG,MyPrint["Integrand LTD1*jac  ="<>ToString[LTDIntegrand[[1]]*jacobian,InputForm]];];
If[DEBUG,MyPrint["Integrand LTD2*jac  ="<>ToString[LTDIntegrand[[2]]*jacobian,InputForm]];];
If[DEBUG,MyPrint["Integrand LTD3*jac  ="<>ToString[LTDIntegrand[[3]]*jacobian,InputForm]];];
If[DEBUG,MyPrint["Integrand UVCT*jac  ="<>ToString[LTDIntegrand[[4]]*jacobian,InputForm]];];
If[DEBUG,MyPrint["Total[LTD]*jac      ="<>ToString[Total[LTDIntegrand]*jacobian,InputForm]];];
Integrand=Integrand*h[rescaling,σ]*((rescaling^6));
If[DEBUG,MyPrint["h[rescaling,σ]="<>ToString[h[rescaling,σ]]];];
If[DEBUG,MyPrint["(rescaling^6)="<>ToString[(rescaling^6)]];];
Integrand=Integrand*(1/rescalingJacobian);
If[DEBUG,MyPrint["1/(rescaling jacobian)="<>ToString[1/rescalingJacobian]];];
Integrand=Integrand*jacobian;
If[DEBUG,MyPrint["jacobian="<>ToString[jacobian]];];
If[DEBUG,MyPrint["Integrand="<>ToString[Integrand]];];
If[DEBUG,MyPrint["Observable argument #1 kmu="<>ToString[kmu,InputForm]];];
If[DEBUG,MyPrint["Observable argument #2 QMinuskmu="<>ToString[QMinuskmu,InputForm]];];
Re[Integrand]
]

BxLIntegrandFunction[{0.1,0.2,0.3,0.999,0.5,0.6},DEBUG->True]

Real-emission A (RxRA)

Polynomial Numerator

RxRANumeratorPolynomial=((-16*I)*ge^4*gs^2*(2*(kx*lx+ky*ly)*(-d^2-2*(8+l0+kx*lx+ky*ly)+d*(8+l0+kx*lx+ky*ly))-2*(-2+d)*kz*(-6+d-kx*lx-ky*ly)*lz+(-2+d)*k0^2*(-2*l0+l0^2-lx^2-ly^2-lz^2)-(-2+d)*kz^2*(-2*l0+l0^2-lx^2-ly^2-lz^2)-(-2+d)*k0*((-2+d)*l0^2-2*kx*lx+4*lx^2-d*lx^2-2*ky*ly+4*ly^2-d*ly^2+l0*(4-2*d+2*kx*lx+2*ky*ly)+2*lz^2-d*lz^2))*(-1+Nc^2)*Qd^2*Qe^2*TF)/.NumericalValues

RxRAIntegrand[k_List,l_List]:=Block[{p1=SelectedP1,p2=SelectedP2},GeVtoPbConversion(((1/(LDot[p1+p2,p1+p2]))(1/(LDot[k-p1-p2,k-p1-p2])))(RxRANumeratorPolynomial/.{
k0->k[[1]],kx->k[[2]],ky->k[[3]],kz->k[[4]],
l0->l[[1]],lx->l[[2]],ly->l[[3]],lz->l[[4]]
})((1/(LDot[p1+p2,p1+p2]))(1/(LDot[l,l])))(1/(((2π))^8))(1/(2LDot[p1+p2,p1+p2]))(((-2π I))/(2Abs[Energy[k]]))(((-2π I))/(2Abs[Energy[k]-Energy[l]]))(((-2π I))/(2Abs[(Energy[p1]+Energy[p2])-Energy[l]])))
]

RxRAIntegrand[{3,4,5,6},{7,8,9,10}]

Defining the complete integrand

RxRARescalingSolutionParametric=Solve[
FullSimplify[Sqrt[(t^2)INldotl]+Sqrt[(t^2)INlkdotlk]+Sqrt[(t^2)INkdotk]-INQ0==0,
Assumptions->{t>0,INkdotk>0,INlkdotlk>0,INldotl>0,INQ0>0}],{t}
];
RxRARescalingSolutionParametricJacobian=
Abs[(D[FullSimplify[Sqrt[(t^2)INldotl]+Sqrt[(t^2)INlkdotlk]+Sqrt[(t^2)INkdotk]-INQ0,
Assumptions->{t>0,INkdotk>0,INlkdotlk>0,INldotl>0,INQ0>0}],t])/.{t->RxRARescalingSolutionParametric}];
RxRARescalingSolution[Q0_,kvec_,lvec_]:=(RxRARescalingSolutionParametric/.{INldotl->VDot[-lvec,-lvec],INlkdotlk->VDot[kvec-lvec,kvec-lvec],INkdotk->VDot[kvec,kvec],INQ0->Q0})
RxRARescalingSolutionJacobian[Q0_,kvec_,lvec_]:=(RxRARescalingSolutionParametricJacobian/.{INldotl->VDot[-lvec,-lvec],INlkdotlk->VDot[kvec-lvec,kvec-lvec],INkdotk->VDot[kvec,kvec],INQ0->Q0})

RxRAIntegrandFunction[xs_,OptionsPattern[DEBUG->False]]:=Module[
{
kmap=SphericalMap[xs[[1]],xs[[2]],xs[[3]]],
(*kmap=CartesianMap[xs[[1]],xs[[2]],xs[[3]]],*)
lmap=SphericalMap[xs[[4]],xs[[5]],xs[[6]]],
(*lmap=CartesianMap[xs[[4]],xs[[5]],xs[[6]]],*)
jacobian=1,Integrand=1,kvec,lvec,kmu,lmu,k0,l0,rescaling,rescalingJacobian,Q0,σ,
DEBUG=OptionValue[DEBUG],lmk,Qml,p1=SelectedP1,p2=SelectedP2
},
Q0=Energy[p1]+Energy[p2];
σ=Q0;
If[DEBUG,MyPrint["xs="<>ToString[xs]];];
lvec=lmap["vector"];
jacobian=jacobian*lmap["jacobian"];
kvec=kmap["vector"];
jacobian=jacobian*kmap["jacobian"];
rescaling=(t/.RxRARescalingSolution[Q0,lvec,kvec][[1]]);
If[DEBUG,MyPrint["rescaling="<>ToString[rescaling]];];
rescalingJacobian=RxRARescalingSolutionJacobian[Q0,lvec,kvec];
lvec=rescaling*lvec;
kvec=rescaling*kvec;
k0=Sqrt[VDot[kvec,kvec]];
(* Pick the negative energy solution *)
l0=-Sqrt[VDot[lvec,lvec]]+Q0;
kmu={k0,kvec[[1]],kvec[[2]],kvec[[3]]};
lmu={l0,lvec[[1]],lvec[[2]],lvec[[3]]};
lmk=lmu-kmu;
Qml={Q0-l0,-lvec[[1]],-lvec[[2]],-lvec[[3]]};
If[DEBUG,MyPrint["lmu=("<>ToString[l0]<>","<>ToString[lvec]<>")"];];
If[DEBUG,MyPrint["kmu=("<>ToString[k0]<>","<>ToString[kvec]<>")"];];
If[DEBUG,MyPrint["-lmu+kmu=("<>ToString[-l0+k0]<>","<>ToString[-lvec+kvec]<>")"];];
If[DEBUG,MyPrint["|l|+|l-k|+|k|="<>ToString[Sqrt[VDot[lvec,lvec]]+Sqrt[VDot[lvec-kvec,lvec-kvec]]+Sqrt[VDot[kvec,kvec]]]<>" vs Q0="<>ToString[Q0]]];
If[DEBUG,MyPrint["|l-k|="<>ToString[Sqrt[VDot[lvec-kvec,lvec-kvec]]]<>" vs (l0-k0)="<>ToString[(l0-k0)]]];
If[DEBUG,MyPrint["|k|="<>ToString[Sqrt[VDot[kvec,kvec]]]<>" vs k0="<>ToString[k0]]];
If[DEBUG,MyPrint["|l|="<>ToString[Sqrt[VDot[lvec,lvec]]]<>" vs -(l0-Q0)="<>ToString[-(l0-Q0)]]];
Integrand=Integrand*RxRAIntegrand[kmu,lmu];
If[DEBUG,MyPrint["Integrand="<>ToString[Integrand]];];
Integrand=Integrand*h[rescaling,σ]*((rescaling^6));
If[DEBUG,MyPrint["h[rescaling,σ]="<>ToString[h[rescaling,σ]]];];
If[DEBUG,MyPrint["(rescaling^6)="<>ToString[(rescaling^6)]];];
Integrand=Integrand*(1/rescalingJacobian);
If[DEBUG,MyPrint["1/(rescaling jacobian)="<>ToString[1/rescalingJacobian]];];
Integrand=Integrand*jacobian;
If[DEBUG,MyPrint["jacobian="<>ToString[jacobian]];];
If[DEBUG,MyPrint["Integrand="<>ToString[Integrand]];];
Re[Integrand]
]

Test the integrand now

RxRAIntegrandFunction[{0.1,0.2,0.3,0.4,0.5,0.6},DEBUG->True]

Real-emission B (RxRB)

Polynomial Numerator

RxRBNumeratorPolynomial=((-16*I)*ge^4*gs^2*(8*ky^2*l0-6*d*ky^2*l0+d^2*ky^2*l0+4*kz^2*l0-4*d*kz^2*l0+d^2*kz^2*l0+2*ky^2*l0^2-d*ky^2*l0^2+2*kz^2*l0^2-d*kz^2*l0^2-32*ky*ly+16*d*ky*ly-2*d^2*ky*ly-4*ky*l0*ly+2*d*ky*l0*ly-4*ky^2*ly^2+2*d*ky^2*ly^2-24*kz*lz+16*d*kz*lz-2*d^2*kz*lz-4*ky*kz*ly*lz+2*d*ky*kz*ly*lz-2*ky^2*lz^2+d*ky^2*lz^2-2*kz^2*lz^2+d*kz^2*lz^2+(-2+d)*k0^2*(-((-2+d)*l0)+l0^2-lz^2)+(-2+d)*kx^2*((-4+d)*l0-l0^2+2*lx^2+lz^2)+2*(-2+d)*k0*(-l0^2+kx*lx+ky*ly+l0*(-2+d-kx*lx-ky*ly)+lz^2)-2*kx*lx*(d^2+2*(8+l0+2*ky*ly+kz*lz)-d*(8+l0+2*ky*ly+kz*lz)))*(-1+Nc^2)*Qd^2*Qe^2*TF)/.NumericalValues

RxRBIntegrand[k_List,l_List]:=Block[{p1=SelectedP1,p2=SelectedP2},GeVtoPbConversion(((1/(LDot[p1+p2,p1+p2]))(1/(LDot[l-p1-p2,l-p1-p2])))(RxRBNumeratorPolynomial/.{
k0->k[[1]],kx->k[[2]],ky->k[[3]],kz->k[[4]],
l0->l[[1]],lx->l[[2]],ly->l[[3]],lz->l[[4]]
})((1/(LDot[p1+p2,p1+p2]))(1/(LDot[k,k])))(1/(((2π))^8))(1/(2LDot[p1+p2,p1+p2]))(((-2π I))/(2Abs[Energy[l]]))(((-2π I))/(2Abs[Energy[k]-Energy[l]]))(((-2π I))/(2Abs[(Energy[p1]+Energy[p2])-Energy[k]])))
]

RxRBIntegrand[{3,4,5,6},{7,8,9,10}]

Defining the complete integrand

RxRBRescalingSolutionParametric=Solve[
FullSimplify[Sqrt[(t^2)INldotl]+Sqrt[(t^2)INlkdotlk]+Sqrt[(t^2)INkdotk]-INQ0==0,
Assumptions->{t>0,INkdotk>0,INlkdotlk>0,INldotl>0,INQ0>0}],{t}
];
RxRBRescalingSolutionParametricJacobian=
Abs[(D[FullSimplify[Sqrt[(t^2)INldotl]+Sqrt[(t^2)INlkdotlk]+Sqrt[(t^2)INkdotk]-INQ0,
Assumptions->{t>0,INkdotk>0,INlkdotlk>0,INldotl>0,INQ0>0}],t])/.{t->RxRBRescalingSolutionParametric}];
RxRBRescalingSolution[Q0_,kvec_,lvec_]:=(RxRBRescalingSolutionParametric/.{INldotl->VDot[lvec,lvec],INlkdotlk->VDot[kvec-lvec,kvec-lvec],INkdotk->VDot[-kvec,-kvec],INQ0->Q0})
RxRBRescalingSolutionJacobian[Q0_,kvec_,lvec_]:=(RxRBRescalingSolutionParametricJacobian/.{INldotl->VDot[lvec,lvec],INlkdotlk->VDot[kvec-lvec,kvec-lvec],INkdotk->VDot[-kvec,-kvec],INQ0->Q0})

RxRBIntegrandFunction[xs_,OptionsPattern[DEBUG->False]]:=Module[
{
kmap=SphericalMap[xs[[1]],xs[[2]],xs[[3]]],
(*kmap=CartesianMap[xs[[1]],xs[[2]],xs[[3]]],*)
lmap=SphericalMap[xs[[4]],xs[[5]],xs[[6]]],
(*lmap=CartesianMap[xs[[4]],xs[[5]],xs[[6]]],*)
jacobian=1,Integrand=1,kvec,lvec,kmu,lmu,k0,l0,rescaling,rescalingJacobian,Q0,σ,
DEBUG=OptionValue[DEBUG],kml,Qmk, p1=SelectedP1, p2=SelectedP2
},
Q0=Energy[p1]+Energy[p2];
σ=Q0;
If[DEBUG,MyPrint["xs="<>ToString[xs]];];
lvec=lmap["vector"];
jacobian=jacobian*lmap["jacobian"];
kvec=kmap["vector"];
jacobian=jacobian*kmap["jacobian"];
rescaling=(t/.RxRBRescalingSolution[Q0,lvec,kvec][[1]]);
If[DEBUG,MyPrint["rescaling="<>ToString[rescaling]];];
rescalingJacobian=RxRBRescalingSolutionJacobian[Q0,lvec,kvec];
lvec=rescaling*lvec;
kvec=rescaling*kvec;
l0=Sqrt[VDot[lvec,lvec]];
(*Must pick negative solution*)
k0=-Sqrt[VDot[kvec,kvec]]+Q0;
kmu={k0,kvec[[1]],kvec[[2]],kvec[[3]]};
lmu={l0,lvec[[1]],lvec[[2]],lvec[[3]]};
kml=kmu-lmu;
Qmk={Q0-k0,-kvec[[1]],-kvec[[2]],-kvec[[3]]};
If[DEBUG,MyPrint["lmu=("<>ToString[l0]<>","<>ToString[lvec]<>")"];];
If[DEBUG,MyPrint["kmu=("<>ToString[k0]<>","<>ToString[kvec]<>")"];];
If[DEBUG,MyPrint["-lmu+kmu=("<>ToString[-l0+k0]<>","<>ToString[-lvec+kvec]<>")"];];
If[DEBUG,MyPrint["|l|+|l-k|+|k|="<>ToString[Sqrt[VDot[lvec,lvec]]+Sqrt[VDot[lvec-kvec,lvec-kvec]]+Sqrt[VDot[kvec,kvec]]]<>" vs Q0="<>ToString[Q0]]];
If[DEBUG,MyPrint["|l-k|="<>ToString[Sqrt[VDot[lvec-kvec,lvec-kvec]]]<>" vs -(l0-k0)="<>ToString[-(l0-k0)]]];
If[DEBUG,MyPrint["|l|="<>ToString[Sqrt[VDot[lvec,lvec]]]<>" vs l0="<>ToString[l0]]];
If[DEBUG,MyPrint["|k|="<>ToString[Sqrt[VDot[kvec,kvec]]]<>" vs k0="<>ToString[k0]]];
If[DEBUG,MyPrint["|k|="<>ToString[Sqrt[VDot[kvec,kvec]]]<>" vs -(k0-Q0)="<>ToString[-(k0-Q0)]]];
Integrand=Integrand*RxRBIntegrand[kmu,lmu];
If[DEBUG,MyPrint["Integrand="<>ToString[Integrand]];];
Integrand=Integrand*h[rescaling,σ]*((rescaling^6));
If[DEBUG,MyPrint["h[rescaling,σ]="<>ToString[h[rescaling,σ]]];];
If[DEBUG,MyPrint["(rescaling^6)="<>ToString[(rescaling^6)]];];
Integrand=Integrand*(1/rescalingJacobian);
If[DEBUG,MyPrint["1/(rescaling jacobian)="<>ToString[1/rescalingJacobian]];];
Integrand=Integrand*jacobian;
If[DEBUG,MyPrint["jacobian="<>ToString[jacobian]];];
If[DEBUG,MyPrint["Integrand="<>ToString[Integrand]];];
Re[Integrand]
]

Test the integrand now

RxRBIntegrandFunction[{0.1,0.2,0.3,0.4,0.5,0.6},DEBUG->True]

Numerical integration of the Double-Triangle supergraph

The integration with Mathematica is not really great, but these expressions do indeed reproduce the expected result around 1011 [pb] when pumping up the stats in Rust.

ThisIntegrand[x1,x2,x3,x4,x5,x6]=LxBIntegrandFunction[{x1,x2,x3,x4,x5,x6},DEBUG->False]+
BxLIntegrandFunction[{x1,x2,x3,x4,x5,x6},DEBUG->False]+
RxRAIntegrandFunction[{x1,x2,x3,x4,x5,x6},DEBUG->False]+
RxRBIntegrandFunction[{x1,x2,x3,x4,x5,x6},DEBUG->False];
(*MCPoints={};*)
rBoundMin=0.01;rBoundMax=0.99;
BoundMin=0.0;BoundMax=1.0;
IntermediateResult=0;
NLOTriangleResults=Association[Table[
IntermediateResult=NIntegrate[
ThisIntegrand[x1,x2,x3,x4,x5,x6],
{x1,rBoundMin,rBoundMax},{x2,BoundMin,BoundMax},{x3,BoundMin,BoundMax},
{x4,rBoundMin,rBoundMax},{x5,BoundMin,BoundMax},{x6,BoundMin,BoundMax},
Method->"AdaptiveMonteCarlo",PrecisionGoal->5,AccuracyGoal->5,MaxPoints->statistics,MaxRecursion->10000
(*,EvaluationMonitor:>AppendTo[MCPoints,{{x1,x2,x3,x4,x5,x6},RealCutRescalingIntegrand[{x1,x2,x3,x4,x5,x6},0,2,False],VirtualCutRescalingIntegrand[{x1,x2,x3,x4,x5,x6},0,2,False]}]*)
];
Print["Result for "<>ToString[statistics]<>" sample points: "<>ToString[IntermediateResult]];
statistics->IntermediateResult
,{statistics, {10000,100000,500000,1000000}}]
]

(*NLOTriangleResults=<|10000->845.4004303717252,100000->879.4715339400847,500000->925.4408263870366,1000000->942.4793227334869,3000000->956.6029873462251,5000000->963.0694576240127,7000000->965.4487212248251,10000000->974.3858074871224,50000000->984.585542008315,100000000->985.8452258234571,500000000->993.8989860249346,1000000000->994.6674869973305|>;*)

FinalNLOTriangleResult= NLOTriangleResults[Max[Keys[NLOTriangleResults]]]

Cross-section in Rust: 1011.14 +- 0.9948 [pb]

NLOTriangleResultBenchmark=1011.0;

NLOVBubbleResultRelDiff=Abs[(FinalNLOTriangleResult-NLOTriangleResultBenchmark)/NLOTriangleResultBenchmark]*100.0

Combining the results into ((Subscript[α, s])/π)!

Compute the relative difference of the total NLO correction w.r.t to the target which is ((Subscript[α, s])/π)Born

Target result:

targetNLOResult=(((((((gs^2)/(4π))))/π))LOResultBenchmark)/.NumericalValues

NLO contribution from the triangle super graph:

FinalNLOTriangleResult

NLO contribution from one bubble graph (there are two):

FinalNLOVBubbleResult

Absolute difference between Final NLO result and target:

(FinalNLOTriangleResult+2*FinalNLOVBubbleResult)-targetNLOResult

Relative difference between Final NLO result and target:

MyPrint["Relative difference: "<>ToString[Abs[((((FinalNLOTriangleResult+2*FinalNLOVBubbleResult)-targetNLOResult)/targetNLOResult))]*100.,InputForm]<>" %."];

The large difference above comes from the very low statistics employed here. If instead we take the rust results (labelled as “benchmark”) we get:

MyPrint["Relative difference: "<>ToString[Abs[((((NLOTriangleResultBenchmark+2*NLOVBubbleResultBenchmark)-targetNLOResult)/targetNLOResult))]*100.,InputForm]<>" %."];

One can also of course directly integrate all contributions at once:

ThisIntegrand[x1,x2,x3,x4,x5,x6]=LxBIntegrandFunction[{x1,x2,x3,x4,x5,x6},DEBUG->False]+
BxLIntegrandFunction[{x1,x2,x3,x4,x5,x6},DEBUG->False]+
RxRAIntegrandFunction[{x1,x2,x3,x4,x5,x6},DEBUG->False]+
RxRBIntegrandFunction[{x1,x2,x3,x4,x5,x6},DEBUG->False]+2(
RBubbleIntegrandFunction[{x1,x2,x3,x4,x5,x6},DEBUG->False]+
alaLTDSoperLRDifferentiatedVBubbleIntegrandFunction[{x1,x2,x3,x4,x5,x6},DEBUG->False]);
(*MCPoints={};*)
rBoundMin=0.001;rBoundMax=0.999;
BoundMin=0.0;BoundMax=1.0;
IntermediateResult=0;
AllNLOResults=Association[Table[
IntermediateResult=NIntegrate[
ThisIntegrand[x1,x2,x3,x4,x5,x6],
{x1,rBoundMin,rBoundMax},{x2,BoundMin,BoundMax},{x3,BoundMin,BoundMax},
{x4,rBoundMin,rBoundMax},{x5,BoundMin,BoundMax},{x6,BoundMin,BoundMax},
Method->"AdaptiveMonteCarlo",PrecisionGoal->5,AccuracyGoal->5,MaxPoints->statistics,MaxRecursion->10000
(*,EvaluationMonitor:>AppendTo[MCPoints,{{x1,x2,x3,x4,x5,x6},RealCutRescalingIntegrand[{x1,x2,x3,x4,x5,x6},0,2,False],VirtualCutRescalingIntegrand[{x1,x2,x3,x4,x5,x6},0,2,False]}]*)
];
Print["Result for "<>ToString[statistics]<>" sample points: "<>ToString[IntermediateResult]];
statistics->IntermediateResult
,{statistics, { 10000,100000,500000}}]
]

Below are hard-coded result that are representative of the integrals you obtained with higher statistics. The Rust implementation fares *much* better of course.

AllNLOResults=<|10000->266.3279180198415,100000->189.80010631798697,500000->224.46257769782477,1000000->232.26211652496414,3000000->251.6080209740919,5000000->255.6637907734908,7000000->258.4061992689343,10000000->266.6870103226783,50000000->273.8477491890012,100000000->276.8124408548081,500000000->283.16369289028125,1000000000->284.6220710686282|>;

FinalNLOResult = AllNLOResults[Max[Keys[AllNLOResults]]]

MyPrint["Relative difference: "<>ToString[Abs[(((FinalNLOResult-targetNLOResult)/targetNLOResult))]*100.,InputForm]<>" %."];