(* ::Package:: *)

(* ::Title:: *)
(*Exercise 3.5*)


(* ::Subtitle:: *)
(*Reducibility of the pentagon in four dimensions*)


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


(* ::Section:: *)
(*Kinematics*)


(* ::Text:: *)
(*We begin by parametrising the kinematics.*)
(*We shall use the symbol "dot[pi,pj]" for the scalar product of two momenta pi and pj. It is convenient to make it commutative as*)


SetAttributes[dot,Orderless];


(* ::Text:: *)
(*We define the external momenta, the on-shell conditions and momentum-conservation constraints on the scalar products, and the chosen set of independent kinematic invariants {s12,s23,s34,s45,s51}. *)


ExternalMomenta = {p1,p2,p3,p4,p5};

OnshellConditions = Table[dot[pi,pi]==0,{pi,ExternalMomenta}];

MomentumConservation = Flatten@Table[Sum[dot[pi,pj],{pi,ExternalMomenta}]==0,{pj,ExternalMomenta}];

InvariantsDefinition = {s12==2*dot[p1,p2], s23==2*dot[p2,p3], s34==2*dot[p3,p4], s45==2*dot[p4,p5], s51==2*dot[p5,p1]};


(* ::Text:: *)
(*We put all the constraints above together and solve them to express all the scalar products of external momenta in terms of the independent kinematic invariants.*)


Kinematics = Flatten@Solve[Join[OnshellConditions,MomentumConservation,InvariantsDefinition]];
Print[Kinematics];


(* ::Section:: *)
(*Gram determinant*)


(* ::Text:: *)
(*We define a function which computes the determinant of the Gram matrix for a given set of vectors.*)


GramDet[q_List]:=Det[Table[dot[qi,qj],{qi,q},{qj,q}]];


(* ::Text:: *)
(*We can use the parametrisation of the kinematics in terms of independent invariants for a self-consistency check: the Gram determinant of {p1,p2,p3,p4,p5} must vanish, as the external momenta are related by momentum conservation. *)


Print["det G(p1,p2,p3,p4,p5) = ", 
  Simplify[GramDet[ExternalMomenta]/.Kinematics]];


(* ::Text:: *)
(*We define the Gram determinant which we will use to reduce the pentagon integral.*)


gramdet=GramDet[{k,p1,p2,p3,p4}];


(* ::Text:: *)
(*As a further self-consistency check, we can verify that this Gram determinant vanishes for a four-dimensional parametrisation of the loop momentum. For instance, we may expand k in {p1,p2,p3,p4}. For this purpose, we need to define a few more rules to implement the bi-linearity of the scalar product.*)


k4D = a1*p1+a2*p2+a3*p3+a4*p4;

ScalarProductRules = {
  dot[a_,b_+c_]:>dot[a,b]+dot[a,c],
  dot[c_*a_,b_]:>c*dot[a,b]/;MemberQ[{a1,a2,a3,a4},c]||NumberQ[c]
  };

Print["det G(k,p1,p2,p3,p4)|(D=4) = ", 
  Simplify[gramdet /. k->k4D //. ScalarProductRules /. Kinematics]];


(* ::Section:: *)
(*Reduction of the pentagon*)


(* ::Text:: *)
(*We compute the rules to express the scalar products involving the loop momentum in terms of the inverse propagators of the pentagon. We denote the latter by d[i] with i=1,...,5.*)


LoopMomenta = {k,k-p1,k-p1-p2,k-p1-p2-p3,k-p1-p2-p3-p4};
InversePropagators = Simplify[Map[dot[#,#]&,LoopMomenta] //.ScalarProductRules/.Kinematics];

LoopScalarProducts = {dot[k,k],dot[k,p1],dot[k,p2],dot[k,p3],dot[k,p4]};

LoopScalarProductsToProp = Flatten@Solve[Thread[Array[d,5]==InversePropagators],LoopScalarProducts];
Print[LoopScalarProductsToProp];


(* ::Text:: *)
(*We can now rewrite the Gram-determinant constraint in terms of inverse propagators and independent kinematic invariants.*)


gramdetInversePropagators = Simplify[gramdet/.LoopScalarProductsToProp/.Kinematics];


(* ::Text:: *)
(*With some manipulations, we can rewrite it in the form 1 + ... = 0.*)


gramdetInversePropagators=Collect[gramdetInversePropagators/(gramdetInversePropagators/.d[_]->0),_d,Simplify];

Print[gramdetInversePropagators];


(* ::Text:: *)
(*In order to write the reduction rule for the pentagon integral explicitly we need to anticipate some notation from chapter 4.*)
(**)
(*We denote by G[a1,a2,a3,a4,a5] the pentagon integral with propagators raised to powers a1,a2,a3,a4,a5. Explicitly, we have*)
(**)
(*G[a1,a2,a3,a4,a5] = \int (-D1)^(-a1) (-D2)^(-a2) (-D3)^(-a3) (-D4)^(-a4) (-D5)^(-a5) .*)
(**)
(*Pay attention to the minus signs in front of the Di's, which are due to our conventions for the Feynman integrals. *)
(*The scalar pentagon integral F_5 in this notation is G[1,1,1,1,1]. G[1,1,1,1,0], G[1,1,0,1,1]... are box integrals, G[1,1,1,0,0], G[1,1,0,1,0]... are triangle integrals, and so on. A negative propagator power is equivalent to a numerator. For instance, G[1,1,1,1,-1] is a box integral G[1,1,1,1,0] with numerator -(k-p1-p2-p3-p4)^2 . *)
(**)
(*Armed with this notation, we multiply the relation above by G[1,1,1,1,1] at the integrand level, and apply the replacement rules of the form*)
(**)
(*G[1,1,1,1,1]*d[5] -> - G[1,1,1,1,0],*)
(*G[1,1,1,1,1]*d[1]*d[3] -> + G[0,1,0,1,1],*)
(**)
(*et cetera, to absorb all the inverse propagators in the definition of the integrals.*)


relation = Expand[gramdetInversePropagators*G[1,1,1,1,1]] //. d[i_]^k_.*G[a__]:>(-1)^k*G@@({a}-k*IdentityMatrix[5][[i]]);

PentagonReduction = G[1,1,1,1,1]/.Flatten[Solve[relation==0,G[1,1,1,1,1]]];
PentagonReduction = Collect[PentagonReduction,_G,Simplify];

Print["F_5 = ", PentagonReduction];
