(* ::Package:: *)

(* ::Title:: *)
(*Exercise 4.16*)


(* ::Subtitle:: *)
(*The massless one-loop box integrals with the differential equations method*)


(* ::Text:: *)
(*This notebook was tested with Wolfram Mathematica 12.3.1.0 and LiteRed v1.83.*)
(*Get LiteRed at https://www.inp.nsk.su/~lee/programs/LiteRed/ .*)


(* ::Text:: *)
(*We start off by setting the working directory and loading LiteRed.*)


$Directory = If[$FrontEnd=!=Null,NotebookDirectory[],DirectoryName[$InputFileName]];
SetDirectory[$Directory];


<<LiteRed`;


(* ::Text:: *)
(*If you have issues loading LiteRed, make sure that it is visible in Mathematica's search path ($Path). If not, add the following line to the file $UserBaseDirectory/Kernel/init.m:*)
(**)
(*$Path = Append[$Path, directory];*)
(**)
(*where the string directory specifies the folder where the file "LiteRed.m" is located. Then quit and restart.*)


(* ::Chapter:: *)
(*a) Defining the integral family and finding a basis of master integrals*)


(* ::Text:: *)
(*Set the spacetime dimension D to 4-2*eps.*)


SetDim[4-2*eps];


(* ::Text:: *)
(*Declare the symbols which represent vectors and kinematic variables.*)


Declare[{k,p1,p2,p3,p4},Vector];
Declare[{s,t},Number];


(* ::Text:: *)
(*Parameterise the kinematics by expressing all scalar products of external momenta "sp[p,q]" in terms of s and t. *)
(*We do not need to consider the scalar products involving p4 because we eliminate the latter through momentum conservation (p4=-p1-p2-p3).*)


sp[p1,p1]=0;
sp[p2,p2]=0;
sp[p3,p3]=0;
sp[p1,p2]=s/2;
sp[p2,p3]=t/2;
sp[p1,p3]=-(s+t)/2;


(* ::Text:: *)
(*Define the integral family with the function NewBasis. *)
(*The first argument is a label to denote the family. *)
(*The second argument is the list of propagators (sp[p] is a shortcut for sp[p,p], and the i0 prescription is absorbed into s and t). *)
(*The third argument is the list of loop momenta. *)
(*The option Directory specifies a directory to store data.*)


NewBasis[
  box,
  {-sp[k], -sp[k+p1], -sp[k+p1+p2], -sp[k+p1+p2+p3]},
  {k},
  Directory->"cache"
];


(* ::Text:: *)
(*Generate the IBP relations.*)


GenerateIBP[box];


(* ::Text:: *)
(*Determine the list of sectors, and which sectors vanish.*)


AnalyzeSectors[box];


(* ::Text:: *)
(*The zero sectors correspond to scaleless integrals, which vanish in dimensional regularisation.*)


Print[ZeroSectors[box]];


(* ::Text:: *)
(*There are 7 non-zero sectors: 2 bubbles (in the s and the t channel), 4 triangles, and 1 box.*)


Print[NonZeroSectors[box]];


(* ::Text:: *)
(*Some sectors may be related. The following command determines which sectors are unique and which are mapped.*)


FindSymmetries[box];


(* ::Text:: *)
(*In this case all sectors are unique.*)


Print[UniqueSectors[box]===NonZeroSectors[box]];


(* ::Text:: *)
(*Solve the IBP relations to reduce the integrals in each of the unique sectors.*)


Map[SolvejSector, UniqueSectors[box]];


(* ::Text:: *)
(*There are 3 master integrals: t-channel bubble, s-channel bubble, and box.*)


g = MIs[box];
nMIs = Length[g];
Print[g];


(* ::Text:: *)
(*By manually attaching a graph to the top sector of the family, we can use LiteRed to draw the integrals.*)


AttachGraph[js[box,1,1,1,1],{
{4->1,"1"},
{1->2,"2"},
{2->3,"3"},
{3->4,"4"},
{1->0,"p1"},
{2->0,"p2"},
{3->0,"p3"},
{4->0,"p4"}
}];


(* ::Text:: *)
(*The graphs of the master integrals are*)


Print[Map[GraphPlot[jGraph[#]]&,g]];


(* ::Text:: *)
(*The IBP reduction of an arbitrary integral of the family is obtained as follows:*)


Print[IBPReduce[ j[box,-1,2,0,1] ]];


(* ::Chapter:: *)
(*b) Differential equations*)


(* ::Text:: *)
(*Compute the derivatives of the master integrals g w.r.t. s and t.*)


gder[s] = Dinv[g,s];
gder[t] = Dinv[g,t];


(* ::Text:: *)
(*IBP-reduce the derivatives to master integrals.*)


gder[s] = IBPReduce[gder[s]];
gder[t] = IBPReduce[gder[t]];


(* ::Text:: *)
(*Extract the matrices of the DEs in the form*)
(**)
(*\[PartialD]g/\[PartialD]s = A[s] . g ,*)
(*\[PartialD]g/\[PartialD]t = A[t] . g .*)


MatrixA[s]=Table[Coefficient[gder[s][[i]],g],{i,nMIs}];
MatrixA[t]=Table[Coefficient[gder[t][[i]],g],{i,nMIs}];

Print["A[s] = ", MatrixForm[MatrixA[s]]];
Print["A[t] = ", MatrixForm[MatrixA[t]]];


(* ::Text:: *)
(*Check the integrability conditions and the scaling relation.*)


IntegrabilityConditions = Simplify[MatrixA[s] . MatrixA[t]-MatrixA[t] . MatrixA[s]+D[MatrixA[s],t]-D[MatrixA[t],s]];

Print["integrability conditions: ",
  DeleteDuplicates[Flatten[IntegrabilityConditions]]==={0}];
  
Scaling = Simplify[s*MatrixA[s]+t*MatrixA[t]];
Print["scaling: ", MatrixForm[Scaling]];

Clear[Scaling, IntegrabilityConditions];


(* ::Chapter:: *)
(*c) Canonical differential equations*)


(* ::Text:: *)
(*Define the basis of pure master integrals f.*)


f = c*{
  s*t*j[box,1,1,1,1],
  s*j[box,1,1,1,0],
  t*j[box,1,1,0,1]
  };


(* ::Text:: *)
(*The overall kinematic-independent factor c is irrelevant for the DEs. We will need its value in point e, for computing the boundary values.*)


PrefactorRule = c->Exp[eps*EulerGamma]*eps^2;


(* ::Text:: *)
(*In order to compute the transformation matrix T from f to g (g = T.f), we IBP reduce f to g.*)


fReduced = IBPReduce[f];

T=Factor[Inverse[Table[Coefficient[mi,g],{mi,fReduced}]]];

Print["T^(-1) = ", MatrixForm[Inverse[T]]];


(* ::Text:: *)
(*Compute the derivatives of f w.r.t. s and t, and IBP reduce.*)


fder[s] = IBPReduce[Dinv[f, s]];
fder[t] = IBPReduce[Dinv[f, t]];


(* ::Text:: *)
(*The derivatives of f are written in terms of the integral basis g. In other words, we have*)
(**)
(*\[PartialD]f/\[PartialD]s = B'[s].g ,*)
(*\[PartialD]f/\[PartialD]t = B'[t].g .*)


MatrixBprime[s] = Table[Coefficient[fder[s][[i]],g],{i,nMIs}];
MatrixBprime[t] = Table[Coefficient[fder[t][[i]],g],{i,nMIs}];


(* ::Text:: *)
(*Use the transformation matrix T to change basis on the right-hand sides of the above DEs. We obtain*)
(**)
(*\[PartialD]f/\[PartialD]s = eps * B[s].f ,*)
(**)
(*where*)
(**)
(*B[s] = B'[s] . T / eps ,*)
(**)
(*and similarly for t.*)


MatrixB[s] = Factor[MatrixBprime[s] . T/eps];
MatrixB[t] = Factor[MatrixBprime[t] . T/eps];

Print["B[s] = ", MatrixForm[MatrixB[s]]];
Print["B[t] = ", MatrixForm[MatrixB[t]]];


(* ::Text:: *)
(*B[s] and B[t] no longer depend on eps.*)


Print["factorisation of eps: ",
  FreeQ[Simplify[{MatrixB[s],MatrixB[t]}],eps]];


(* ::Text:: *)
(*Alternatively, we may obtain B[s] and B[t] also through the "gauge" transformation*)
(**)
(*eps B[s] = T^(-1) . A[s] . T - T^(-1) . \[PartialD]T/\[PartialD]s ,*)
(**)
(*and similarly for t.*)


Print["check gauge transformation: ",
  And[
    Simplify[Inverse[T] . MatrixA[s] . T-Inverse[T] . D[T,s]] === Simplify[eps*MatrixB[s]],
    Simplify[Inverse[T] . MatrixA[t] . T-Inverse[T] . D[T,t]] === Simplify[eps*MatrixB[t]]
  ]];


(* ::Text:: *)
(*Check integrability conditions and scaling relation.*)


IntegrabilityConditions = Simplify[{
  MatrixB[s] . MatrixB[t]-MatrixB[t] . MatrixB[s],
  D[MatrixB[s],t]-D[MatrixB[t],s]
  }];
Print[DeleteDuplicates[Flatten[IntegrabilityConditions]]==={0}];

Scaling = Simplify[s*MatrixB[s]+t*MatrixB[t]];
Print["scaling: ", MatrixForm[Scaling]];

Clear[Scaling,IntegrabilityConditions];


(* ::Chapter:: *)
(*d) Change of variables*)


(* ::Text:: *)
(*By changing variables from (s,t) to (s,x) with x = t/s we obtain the following system of DEs,*)
(**)
(*\[PartialD]f/\[PartialD]s = C[s] . f  ,*)
(*\[PartialD]f/\[PartialD]x = C[x] . f ,*)
(**)
(*where*)
(**)
(*C[s] = B[s] + t/s*B[t] ,*)
(*C[x] = s B[t] .*)


MatrixC[s] = Simplify[MatrixB[s]+t/s*MatrixB[t] /. t->s*x];
Print["C[s] = ", MatrixForm[MatrixC[s]]];
    
MatrixC[x] = Simplify[s*MatrixB[t]/. t->s*x];
Print["C[x] = ", MatrixForm[MatrixC[x]]];


(* ::Text:: *)
(*We perform a partial fraction decomposition of the matrix C[x] to highlight the irreducible factors in the denominators, which give the alphabet letters.*)


MatrixC[x] = Apart[MatrixC[x],x];
Print["C[x] = ", MatrixForm[MatrixC[x]]];


(* ::Text:: *)
(*Integrate C[x] and C[s] to write the DE in differential form*)
(**)
(*d f = eps (d Ctilde) . f . *)
(**)
(*Note that this only works because of the separation of variables: C[x] only depend on x, and C[s] only on s.*)


MatrixCtilde = Expand[Integrate[MatrixC[x],x]+Integrate[MatrixC[s],s]];
Print["Ctilde = ", MatrixForm[MatrixCtilde]];


(* ::Text:: *)
(*Check that D[Ctilde,s] = C[s] and D[Ctilde,x] = C[x].*)


Print["check Ctilde: ",
  And[Simplify[D[MatrixCtilde,s]]===MatrixC[s], Simplify[D[MatrixCtilde,x]]===Simplify[MatrixC[x]]]];


(* ::Text:: *)
(*Extract the constant matrices which multiply the logarithms.*)


C1=Coefficient[MatrixCtilde,Log[s]];
C2=Coefficient[MatrixCtilde,Log[x]];
C3=Coefficient[MatrixCtilde,Log[1+x]];

Print["C_1 = ", MatrixForm[C1]];
Print["C_2 = ", MatrixForm[C2]];
Print["C_3 = ", MatrixForm[C3]];


(* ::Chapter:: *)
(*e) Compute the symbol of the solution*)


(* ::Text:: *)
(*We use the known expression of the bubble integrals to obtain the analytic expression of the master integrals of the bubble type.*)


Bubble[a1_,a2_,d_][p2_]:=Gamma[a1+a2-d/2]*Gamma[d/2-a1]*Gamma[d/2-a2]/Gamma[a1]/Gamma[a2]/Gamma[d-a1-a2]*(-p2)^(d/2-a1-a2);

f2 = fReduced[[2]] /. PrefactorRule /. j[box,1,0,1,0]:>Bubble[1,1,4-2*eps][s] /. t->s*x;
f3 = fReduced[[3]] /. PrefactorRule /. j[box,0,1,0,1]:>Bubble[1,1,4-2*eps][t] /. t->s*x;

Print["f_2 = ", Series[f2,{eps,0,2}]];


(* ::Text:: *)
(*The terms of order eps^0 give the weight-0 boundary values. We leave the value of the box-type master integral as a free parameter b[0], and fix it by requiring the finiteness of the solution at x=-1.*)


BdryValuesW0 = {
  b[0],
  Normal[Series[f2,{eps,0,0}]],
  Normal[Series[f3,{eps,0,0}]]
};

BdryValuesW0constraint = Flatten[Solve[C3 . BdryValuesW0 == 0, b[0]]];
BdryValuesW0 = BdryValuesW0 /. BdryValuesW0constraint;

Print["weight-0 boundary values = ", BdryValuesW0]; 
Clear[BdryValuesW0constraint];


(* ::Text:: *)
(*The following function takes as input *)
(*- the connection matrix of the DEs in differential form (written in terms of logarithms), *)
(*- the weight-0 boundary values, *)
(*- the desired maximum weight, *)
(*and returns the symbol of the solution up to the specified weight.*)


GetSymbol[Atilde_List, bdry_List, maxweight_Integer]:=Module[{SBw,w},

  SBw[0] = bdry*SB[]; (* SB[] = 1, weight-0 symbol *)

  For[w=1,w<=maxweight,w++,
    SBw[w] = Expand[Atilde . SBw[w-1]] /. SB[a___]*Log[b_]:>SB[a,b];
  ]; 

  Return[Sum[SBw[w]*eps^w,{w,0,maxweight}]];
];


(* ::Text:: *)
(*Up to weight 4 we obtain*)


MaxWeight=4;
fSB = GetSymbol[MatrixCtilde, BdryValuesW0, MaxWeight];

Print["S(f[1]) = ", fSB[[1]]+O[eps]^(MaxWeight+1)];
Print["S(f[2]) = ", fSB[[2]]+O[eps]^(MaxWeight+1)];
Print["S(f[3]) = ", fSB[[3]]+O[eps]^(MaxWeight+1)];


(* ::Chapter:: *)
(*f) Analytic solution*)


(* ::Text:: *)
(*We remove the overall factor of (-s)^(-eps), and expand the analytic expressions for the bubble-type master integrals around eps=0 up to the chosen maximum weight. In order to simplify the expressions, we specify the domain of all kinematic variables.*)


MaxWeight = 3;

EuclideanRegion={s<0,t<0,x>0};

h2 = Simplify[Normal[Series[(-s)^eps*f2,{eps,0,MaxWeight}]],EuclideanRegion];
h3 = Simplify[Normal[Series[(-s)^eps*f3,{eps,0,MaxWeight}]],EuclideanRegion];

Print["h2 = ", h2];
Print["h3 = ", h3];


(* ::Text:: *)
(*PolyGamma[2,1] is nothing but -2 Zeta[3] , which can be revealed through the function FunctionExpand.*)


h2 = Collect[h2,eps,FunctionExpand];
h3 = Collect[h3,eps,FunctionExpand];

Print["h2 = ", h2];
Print["h3 = ", h3];


(* ::Text:: *)
(*We prepare the boundary values at x0=1. We denote the unknown weight-w value of the box-type integral as e[w].*)


BdryValuesH[0] = BdryValuesW0;

x0=1;

Do[
  BdryValuesH[w] = {
    e[w],
    Coefficient[h2/.x->x0,eps,w],
    Coefficient[h3/.x->x0,eps,w]
 };,
{w,1,MaxWeight}];


(* ::Text:: *)
(*We integrate the DEs for h order by order in eps, and fix the undetermined boundary value by requiring the finiteness of the solution at xspurious=-1. Assuming that x>1 simplifies the integration. The result may then be analytically continued to x<1.*)


Clear[h$];

xspurious=-1;

h$[0] = BdryValuesH[0];

For[w=1, w<=MaxWeight, w++,
  Print["weight ", w];
  
  constraints = DeleteCases[ Limit[ C3 . h$[w-1] , x->xspurious] , 0];
  constraints = Flatten[Solve[constraints==0,e[w-1]]];
  Print["new constraints on the boundary values: ", constraints];

  BdryValuesH[w-1] = BdryValuesH[w-1] /. constraints;
  h$[w-1] = Expand[h$[w-1] /. constraints];
  
  h$[w] = Integrate[ D[MatrixCtilde,x] . h$[w-1]/.x->t,{t,x0,x},Assumptions->x>1]+BdryValuesH[w];
  h$[w] = FunctionExpand[h$[w],EuclideanRegion];
 
  Print["solution: ", h$[w]]; 
  Print["------------------------------------------"];
  Clear[constraints];
];

constraints = DeleteCases[ Limit[ C3 . h$[MaxWeight] , x->xspurious] , 0];
constraints = Flatten[Solve[constraints==0,e[MaxWeight]]];
Print["new constraints on the boundary values: ", constraints];

BdryValuesH[MaxWeight] = BdryValuesH[MaxWeight] /. constraints;
h$[MaxWeight] = Expand[h$[MaxWeight] /. constraints];
  
hAnalytic = Sum[h$[w]*eps^w,{w,0,MaxWeight}];

Clear[h$, constraints];


Print["h[1] = ", hAnalytic[[1]]+O[eps]^(MaxWeight+1)];
Print["h[2] = ", hAnalytic[[2]]+O[eps]^(MaxWeight+1)];
Print["h[3] = ", hAnalytic[[3]]+O[eps]^(MaxWeight+1)];


(* ::Text:: *)
(*Check that the solution to the DEs equals the known analytic expressions for the bubble-type master integrals.*)


Print["check known solutions: ", 
  Simplify[{hAnalytic[[2]]-h2,hAnalytic[[3]]-h3}]==={0,0}];


(* ::Chapter:: *)
(*g) Comparison against Mellin-Barnes result*)


(* ::Text:: *)
(*Hard-coded expression for the box integral F4 = j[box,1,1,1,1] obtained with the Mellin-Barnes method up to order eps^1 (corresponding to MaxWeight = 3):*)


F4MellinBarnes = Exp[-eps*EulerGamma]*(-s)^(-eps)/s/t*
  (4/eps^2
  -2/eps*Log[x]
  -4*Pi^2/3 
  + eps*(2*PolyLog[3,-1/x]+2*Log[x]*PolyLog[2,-1/x]+Log[x]^3+7*Pi^2/6*Log[x]-
    (Log[x]^2+Pi^2)*Log[1+x]-34/3*Zeta[3])) + O[eps]^2;


(* ::Text:: *)
(*Corresponding result from the DEs (MaxWeight must be at least 3):*)


F4DE = (-s)^(-eps)*(T . hAnalytic)[[3]] /. PrefactorRule;
F4DE = Normal[Series[F4DE,{eps,0,1}]];


(* ::Text:: *)
(*The first three eps-orders match identically.*)


zero = F4DE-F4MellinBarnes;

Print["eps^(-2) -> ", Simplify[Coefficient[zero,eps,-2]]];
Print["eps^(-1) -> ", Simplify[Coefficient[zero,eps,-1]]];
Print["eps^0 -> ", Simplify[Coefficient[zero,eps,0]]];


(* ::Text:: *)
(*At order eps^1 we need some functional identities among the polylogarithms:*)


PolylogIdentities = {
  PolyLog[2,-1/x] -> -PolyLog[2,-x]-Log[x]^2/2-Zeta[2],
  PolyLog[3,-1/x] -> PolyLog[3,-x]+Log[x]^3/3!+Zeta[2]*Log[x]
};

Print["eps^1 -> ", Simplify[Coefficient[zero,eps,1]/.PolylogIdentities]];
