(* ::Package:: *)

(* ::Title:: *)
(*Exercise 3.9*)


(* ::Subtitle:: *)
(*Rank-one triangle reduction with the direct extraction method*)


(* ::Text:: *)
(*This notebook was tested with Wolfram Mathematica 12.3.1.0.*)


(* ::Section:: *)
(*Scalar product*)


(* ::Text:: *)
(*We begin by defining a symbol for the scalar product of two vectors. As we did in exercise 3.5, we use "dot[a,b]" for the scalar product of two vectors "a" and "b". We make "dot" commutative, and define a few rules to simplify the expressions.*)
(*We denote by gam[a,b] the light-like vector "1/2*<a|\gamma^{\mu}|b]" made from the negative-helicity spinor "a" and the positive-helicity spinor "b".*)


SetAttributes[dot,Orderless];


ScalarProductRules = {
(* short-hand notation *)
  dot[a_]:>dot[a,a],
  
(* linearity of the scalar product *)
  dot[a_,b_+c_]:>dot[a,b]+dot[a,c],
  dot[c_*a_,b_]:>c*dot[a,b]/;NumberQ[c],
  
(* spinor contractions *)
  dot[gam[a_,b_],gam[b_,a_]]:>-dot[a,b],
  dot[gam[a_,b_],gam[a_,b_]]:>0
};


(* ::Section:: *)
(*Triangle definition*)


(* ::Text:: *)
(*We define the inverse propagators of the triangle with external momenta P, Q and R, as well as the integrand under investigation. We use dot[a] as a short-hand for dot[a,a].*)


D1 = dot[k]-m1sq;
D2 = dot[k-P]-m2sq;
D3 = dot[k-P-Q]-m3sq;

Integrand = dot[k,Z]/(-D1)/(-D2)/(-D3);


(* ::Text:: *)
(*We add to the rules for the scalar product the definition of the triangle kinematic invariants S, T, and U.*)


ScalarProductRules = Join[ScalarProductRules,
  {dot[P,P]->S, dot[Q,Q]->T, dot[R,R]->U}];


(* ::Text:: *)
(*It is also convenient to define the auxiliary variables hatS and hatT (as in part (a) of this exercise).*)


ExtraDefinitions = {
 hatS->S+m1sq-m2sq,
 hatT->T+m2sq-m3sq+2*dot[Q,P]
};


(* ::Section:: *)
(*Numerators*)


(* ::Text:: *)
(*We define the general (four-dimensional) parametrisation of the numerators for the triangle with momenta P,Q and R,*)


Delta["P|Q|R"] = c[0,"P|Q|R"] +
  c[1,"P|Q|R"]*dot[k,\[Omega]1tri]+c[2,"P|Q|R"]*dot[k,\[Omega]2tri]+
  c[3,"P|Q|R"]*(dot[k,\[Omega]1tri]^2-dot[\[Omega]1tri,\[Omega]1tri]/dot[\[Omega]2tri,\[Omega]2tri]*dot[k,\[Omega]2tri]^2) + 
  c[4,"P|Q|R"]*dot[k,\[Omega]1tri]*dot[k,\[Omega]2tri]+
  c[5,"P|Q|R"]*dot[k,\[Omega]1tri]^3+
  c[6,"P|Q|R"]*dot[k,\[Omega]1tri]^2*dot[k,\[Omega]2tri];


(* ::Text:: *)
(*and for the P-channel bubble*)


Delta["P|QR"]=c[0,"P|QR"] +
  c[1,"P|QR"] dot[k,\[Omega]1bub]+
  c[2,"P|QR"]dot[k,\[Omega]2bub]+
  c[3,"P|QR"]dot[k,\[Omega]3bub]+
  c[4,"P|QR"](dot[k,\[Omega]1bub]^2-dot[\[Omega]1bub,\[Omega]1bub]/dot[\[Omega]3bub,\[Omega]3bub]*dot[k,\[Omega]3bub]^2)+
  c[5,"P|QR"](dot[k,\[Omega]2bub]^2-dot[\[Omega]2bub,\[Omega]2bub]/dot[\[Omega]3bub,\[Omega]3bub]*dot[k,\[Omega]3bub]^2)+
  c[6,"P|QR"] dot[k,\[Omega]1bub]dot[k,\[Omega]2bub]+
  c[7,"P|QR"] dot[k,\[Omega]1bub]dot[k,\[Omega]3bub]+
  c[8,"P|QR"] dot[k,\[Omega]2bub]dot[k,\[Omega]3bub];


(* ::Section:: *)
(*Spurious spaces*)


(* ::Text:: *)
(*We define vectors spanning the spurious spaces for the bubble (\[Omega]1bub, \[Omega]2bub, \[Omega]3bub) and the triangle (\[Omega]1tri, \[Omega]2tri). We do so by defining rules for the scalar product between them and an arbitrary vector b. The definitions follow those in the text (with Pf = P^\[Flat], Px = Overscript[P,\[Hacek]], Qx = Overscript[Q,\[Hacek]]). \[CapitalPhi]bub and \[CapitalPhi]tri are arbitrary factors which ensure that all vectors have the correct phase. As well see, we do not need to know their expression, as they will cancel out in the result.*)


SpuriousSpaces = {
  dot[\[Omega]1bub,b_]:>dot[b,gam[Pf,n]]*\[CapitalPhi]bub+dot[b,gam[n,Pf]]/\[CapitalPhi]bub,
  dot[\[Omega]2bub,b_]:>dot[b,gam[Pf,n]]*\[CapitalPhi]bub-dot[b,gam[n,Pf]]/\[CapitalPhi]bub,
  dot[\[Omega]3bub,b_]:>dot[Pf,b]-S*dot[n,b]/2/dot[P,n],
  
  dot[\[Omega]1tri,b_]:>dot[b,gam[Px,Qx]]*\[CapitalPhi]tri+dot[b,gam[Qx,Px]]/\[CapitalPhi]tri,
  dot[\[Omega]2tri,b_]:>dot[b,gam[Px,Qx]]*\[CapitalPhi]tri-dot[b,gam[Qx,Px]]/\[CapitalPhi]tri
};


(* ::Text:: *)
(*We further define some rules to simplify the scalar products involving the spurious vectors.*)


ScalarProductRules = Join[ScalarProductRules,
  {
(* bubble spurious space *)
  dot[n,n]->0,
  dot[Pf,Pf]->0,
  dot[P,Pf]->S/2,
  dot[n,Pf]->dot[n,P],
  
  dot[n,gam[n,Pf]]->0,
  dot[n,gam[Pf,n]]->0,
  
  dot[Pf,gam[n,Pf]]->0,
  dot[Pf,gam[Pf,n]]->0,  
  
  dot[P,gam[n,Pf]]->0,
  dot[P,gam[Pf,n]]->0,

(* triangle spurious space *)    
  dot[Px,Px]->0,
  dot[Qx,Qx]->0,
  dot[P,gam[Qx,Px]]->0,
  dot[Q,gam[Qx,Px]]->0,
  dot[P,gam[Px,Qx]]->0,
  dot[Q,gam[Px,Qx]]->0  
}];


(* ::Text:: *)
(*We can now check that the spurious-basis vectors are orthogonal to the corresponding physical vectors, and among each other.*)


Print["check orthogonality to physical spaces: ",
  DeleteDuplicates[
    {dot[\[Omega]1bub,P],dot[\[Omega]2bub,P],dot[\[Omega]3bub,P],
     dot[\[Omega]1tri,P],dot[\[Omega]2tri,P],dot[\[Omega]1tri,Q],dot[\[Omega]2tri,Q]} //.
       SpuriousSpaces//.ScalarProductRules]==={0}];

Print["check orthogonality of spurious-basis vectors: ", 
  DeleteDuplicates[
    {dot[\[Omega]1tri,\[Omega]2tri],
     dot[\[Omega]1bub,\[Omega]2bub],dot[\[Omega]1bub,\[Omega]3bub],dot[\[Omega]2bub,\[Omega]3bub]}//.
       SpuriousSpaces/.ScalarProductRules]==={0}];


(* ::Text:: *)
(*As a further example, we can compute the norms of the spurious-basis vectors:*)


(* compute norms of spurious-basis vectors *) 
Print["\[Omega]1bub^2 = ", Simplify[dot[\[Omega]1bub,\[Omega]1bub]//.SpuriousSpaces//.ScalarProductRules]];
Print["\[Omega]2bub^2 = ", Simplify[dot[\[Omega]2bub,\[Omega]2bub]//.SpuriousSpaces//.ScalarProductRules]];
Print["\[Omega]3bub^2 = ", Simplify[dot[\[Omega]3bub,\[Omega]3bub]//.SpuriousSpaces//.ScalarProductRules]];
Print["\[Omega]1tri^2 = ", Simplify[dot[\[Omega]1tri,\[Omega]1tri]//.SpuriousSpaces//.ScalarProductRules]];
Print["\[Omega]2tri^2 = ", Simplify[dot[\[Omega]2tri,\[Omega]2tri]//.SpuriousSpaces//.ScalarProductRules]];


(* ::Text:: *)
(*It is convenient to keep separate the rules for the auxiliary vectors Px and Qx. Requiring that they are light-like gives the two possible values for \[Gamma] (Subscript[\[Gamma],+]and Subscript[\[Gamma],-] in the text).*)


PxQxrules = {dot[Px,b_]:>\[Gamma]/(\[Gamma]^2-S*T)*(\[Gamma]*dot[P,b]-S*dot[Q,b]),
		     dot[Qx,b_]:>\[Gamma]/(\[Gamma]^2-S*T)*(\[Gamma]*dot[Q,b]-T*dot[P,b])};
		     
\[Gamma]pm=Solve[Simplify[{dot[Px,Px],dot[Qx,Qx]}/\[Gamma]^2//.PxQxrules//.ScalarProductRules]==0,\[Gamma]];
Print["{", Subscript["\[Gamma]","-"], ", ", Subscript["\[Gamma]","+"], "} = ", \[Gamma]/.\[Gamma]pm];


(* ::Section:: *)
(*Extraction of the triangle coefficients*)


(* ::Text:: *)
(*We hard code the loop-momentum parametrisation on the triple cut P|Q|R (D1,D2,D3=0) discussed in the text.*)


TriangleLoopParametrisation = 
  dot[k,b_]:>b1*dot[P,b]+b2*dot[Q,b]+(t+\[Gamma]tri/t)/2*dot[\[Omega]1tri,b]+(t-\[Gamma]tri/t)/2*dot[\[Omega]2tri,b] /.
    {b1->(hatS*T-dot[P,Q]*hatT)/2/(S*T-dot[P,Q]^2),
     b2->(S*hatT-hatS*dot[P,Q])/2/(S*T-dot[P,Q]^2)};


(* ::Text:: *)
(*We use \[Gamma]tri here to distinguish this parameter from the \[Gamma] used in the projections Px and Qx. We obtain its value by requiring that D1 = 0.*)


\[Gamma]triSol=Simplify[Flatten[Solve[(D1//.ScalarProductRules//.TriangleLoopParametrisation//.SpuriousSpaces//.ScalarProductRules)==0,\[Gamma]tri]]];


(* ::Text:: *)
(*We can now verify that the three inverse propagators indeed vanish with this loop-momentum parametrisation.*)


Print["{D1,D2,D3} = ", 
  Simplify[{D1,D2,D3}//.ScalarProductRules//.
    TriangleLoopParametrisation//.SpuriousSpaces//.ScalarProductRules/.\[Gamma]triSol/.ExtraDefinitions]];


(* ::Text:: *)
(*The coefficients of the box numerator vanish, as all quadruple cuts are zero.*)
(*The triangle coefficient c[0;"P|Q|R"] can be extracted directly using the Inf operation defined in the text.*)


Inf[expr_,x_]:=Normal[Series[expr,{x,Infinity,0}]];


DeltaTriangle = Delta["P|Q|R"] //. TriangleLoopParametrisation //. SpuriousSpaces //. ScalarProductRules;

Print["(Inf[Delta[P|Q|R],t]_{t^0} + Inf[Delta[P|Q|R],1/t]_{t^0})/2 = ",
  (Coefficient[Inf[DeltaTriangle,t],t,0]+Coefficient[(Inf[DeltaTriangle/.t->1/invt,invt]/.invt->1/t),t,0])/2];


(* ::Text:: *)
(*We are however interested in all coefficients of the triangle numerator. We thus equate the triple cut of the integrand to the general parametrisation of the triangle numerator on the triple cut, and solve the equations order by order in t.*)


TripleCutIntegrand = Simplify[Integrand*D1*D2*D3] //. TriangleLoopParametrisation;
TriangleEquations = CoefficientList[(DeltaTriangle+TripleCutIntegrand)*t^3,t];
TriangleCoefficients = Cases[DeltaTriangle,c[__],Infinity];
TriangleCoefficientsRules = Flatten[Solve[TriangleEquations==0,TriangleCoefficients]];

Print["triangle-numerator coefficients: ", TriangleCoefficientsRules];


(* ::Section:: *)
(*Extraction of bubble coefficients*)


(* ::Text:: *)
(*We hard code the loop-momentum parametrisation on the bubble cut D1,D2=0 discussed in the text, and verify that the inverse propagators D1 and D2 indeed vanish.*)


BubbleLoopParametrisation = 
  dot[k,b_]:>a1*dot[Pf,b]+a2*dot[n,b]+a3*dot[gam[Pf,n],b]*\[CapitalPhi]bub+a4*dot[gam[n,Pf],b]/\[CapitalPhi]bub /.
    {
     a1->y,
     a2->(hatS-S*y)/2/dot[P,n],
     a3->t,
     a4->(y*(hatS-S*y)-m1sq)/2/t/dot[P,n]
    };

Print["{D1,D2} = ", 
  Simplify[{D1,D2}//.ScalarProductRules//.BubbleLoopParametrisation//.ScalarProductRules/.ExtraDefinitions]];


(* ::Text:: *)
(*As an example, we compute the scalar products between the bubble spurious basis vectors and the loop momentum.*)


Print["k\[CenterDot]\[Omega]1bub = ", Simplify[dot[k,\[Omega]1bub] /. SpuriousSpaces //. BubbleLoopParametrisation //. ScalarProductRules]];
Print["k\[CenterDot]\[Omega]2bub = ", Simplify[dot[k,\[Omega]2bub] /. SpuriousSpaces //. BubbleLoopParametrisation //. ScalarProductRules]];
Print["k\[CenterDot]\[Omega]3bub = ", Simplify[dot[k,\[Omega]3bub] /. SpuriousSpaces //. BubbleLoopParametrisation //. ScalarProductRules]];


(* ::Text:: *)
(*We parametrise the general bubble numerator on the bubble cut P|QR.*)


DeltaBubble=Delta["P|QR"]//.SpuriousSpaces//.BubbleLoopParametrisation//.ScalarProductRules;


(* ::Text:: *)
(*We verify that the operator P defined in the text extracts directly the bubble coefficient c[0;"P|QR"].*)


Poperator[func_]:=
  Coefficient[Coefficient[func,t,0],y,0]+
    hatS/2/S*Coefficient[Coefficient[func,t,0],y,1]+
      1/3*(hatS^2/S^2-m1sq/S)*Coefficient[Coefficient[func,t,0],y,2];
      
Print["P[Inf[Inf[\[CapitalDelta][P|QR],t],y]] = ", Simplify[Poperator[Inf[Inf[DeltaBubble,t],y]]]];


(* ::Text:: *)
(*We parametrise the integrand on the bubble cut, and apply the operation for the direct extraction to it.*)


BubbleCutIntegrand = Simplify[Simplify[Integrand*D1*D2] //. ScalarProductRules //. BubbleLoopParametrisation //. ScalarProductRules];

Print["P[Inf[Inf[I*D1*D2,t],y]] = ", Poperator[Inf[Inf[BubbleCutIntegrand,t],y]]];


(* ::Text:: *)
(*We parametrise the triangle subtraction term on the double cut P|Q, and apply the direct-extraction operation to it.*)


TriangleSubtraction=Together[Delta["P|Q|R"]/(-D3) /. TriangleCoefficientsRules //. ScalarProductRules //. SpuriousSpaces //. 
  BubbleLoopParametrisation //.ScalarProductRules];
  
PTriangleSubtraction = Poperator[Inf[Inf[TriangleSubtraction,t],y]];

Print["P[Inf[Inf[TriangleSubtraction,t],y]] = ", PTriangleSubtraction];


(* ::Text:: *)
(*We simplify the expression of the triangle subtraction term with some manual replacements. Here, \[Gamma]5 is the fifth Dirac matrix. Thanks to the anti-symmetry of the trace with four gamma matrices and \[Gamma]5, the traces involving \[Gamma]5 drop out. We implement this by sorting the arguments of these traces according to Mathematica's canonical order, and keeping track of the sign with the function Signature. Finally, we apply the rule for the trace of four Dirac matrices, and simplify.*)


PTriangleSubtraction = Simplify[PTriangleSubtraction /. 
  dot[A_,gam[Qx,Px]]*dot[B_,gam[Px,Qx]]:>trm[A,Px,B,Qx]/4 /.
  trm[x__]:>1/2*tr[x]-1/2*tr[x,\[Gamma]5] /. 
  tr[x__,\[Gamma]5]:>Signature[{x}]*tr5@@Append[Sort[{x}],\[Gamma]5] /.
  tr[x1_,x2_,x3_,x4_]:>4*(dot[x1,x2]*dot[x3,x4]-dot[x1,x3]*dot[x2,x4]+dot[x1,x4]*dot[x2,x3])];
  
PTriangleSubtraction = FullSimplify[PTriangleSubtraction//.PxQxrules//.ScalarProductRules];
Print["P[Inf[Inf[TriangleSubtraction,t],y]] = ", PTriangleSubtraction];


(* ::Text:: *)
(*We take the average of the two projections SubPlus[\[Gamma]] and SubMinus[\[Gamma]], as discussed in the text. For this purpose, we need to apply the rules for the Px and Qx vectors.*)


PTriangleSubtractionAverage = FullSimplify[Sum[PTriangleSubtraction //. PxQxrules /.ScalarProductRules /. \[Gamma]pm[[i]] ,{i,1,2}]/2];
Print["1/2 \sum_{", Subscript["\[Gamma]","+"], ", ", Subscript["\[Gamma]","-"], "} P[Inf[Inf[TriangleSubtraction,t],y]] = ", 
  PTriangleSubtractionAverage];


(* ::Text:: *)
(*Finally, we subtract the triangle term from the integrand on the bubble cut. This completes the direct extraction of the triangle bubble coefficient c[0,"P|QR"].*)


c[0,"P|QR"] = Simplify[Poperator[Inf[Inf[BubbleCutIntegrand,t],y]] - PTriangleSubtractionAverage];
Print["c[0,P|QR] = ", c[0,"P|QR"]];
