(* ::Package:: *)

(* ::Title:: *)
(*Exercise 4.12*)


(* ::Subtitle:: *)
(*The massive bubble 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 2-2*eps.*)


SetDim[2-2*eps];


(* ::Text:: *)
(*Declare the symbols which represent vectors,*)


Declare[{k,p},Vector];


(* ::Text:: *)
(*and the symbols representing the mass (squared) in the internal propagators and s*)


Declare[{s,msq},Number];


(* ::Text:: *)
(*Parametrise the kinematics by expressing all scalar products of external momenta in terms of a minimal set of variables. In this case there is only one scalar product:*)


sp[p,p]=s;


(* ::Text:: *)
(*Define the integral family. *)
(*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[
  bubble,
  {-sp[k]+msq, -sp[k+p]+msq},
  {k},
  Directory->"cache"
];


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


GenerateIBP[bubble];


(* ::Text:: *)
(*There are two IBP relations. For generic propagators exponents a1 and a2 they are*)


IBPrelations = Collect[IBP[bubble]@@{a1,a2},_j];
Print[TableForm[IBPrelations]];


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


AnalyzeSectors[bubble];


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


Print[ZeroSectors[bubble]];


(* ::Text:: *)
(*There are 3 non-zero sectors: 2 tadpoles and the bubble.*)


Print[NonZeroSectors[bubble]];


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


FindSymmetries[bubble];


(* ::Text:: *)
(*The tadpole sectors are related. The rules which map one sector to the other are obtained as follows. For example, j[bubble,1,0] is mapped to j[bubble,0,1].*)


Print["mapped sectors: ", MappedSectors[bubble]];
mappings = Flatten[MappedSectors[bubble]/.js->jRules];
Print["rules: ", mappings];
Print["example: ", j[bubble,1,0]->(j[bubble,1,0]/.mappings)];


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


Map[SolvejSector, UniqueSectors[bubble]];


(* ::Text:: *)
(*There are 2 master integrals: a tadpole and a bubble.*)


g = MIs[bubble];
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[bubble,1,1],{
  {2->1,"1"},
  {1->2,"2"},
  {1->0,"p"},
  {2->0,"-p"}
}];


(* ::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[j[bubble,2,0]->IBPReduce[ j[bubble,2,0] ]];
Print[j[bubble,2,1]->IBPReduce[ j[bubble,2,1] ]];


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


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


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

Print["D[g,msq] = ", gder[msq]]


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


gder[s] = IBPReduce[gder[s]];
gder[msq] = IBPReduce[gder[msq]];


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


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

Print["A[s] = ", MatrixForm[Coefficient[MatrixA[s],eps,0]] + eps*MatrixForm[Coefficient[MatrixA[s],eps,1]]];
Print["A[msq] = ", MatrixForm[Coefficient[MatrixA[msq],eps,0]] + eps*MatrixForm[Coefficient[MatrixA[msq],eps,1]]];


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


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

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

Clear[Scaling, IntegrabilityConditions];


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


(* ::Text:: *)
(*Change integral basis from g to f = T . g, with*)


T = {{1,0},{0,Sqrt[(-s)*(4*msq-s)]}};
f = T . g;


(* ::Text:: *)
(*Verify that the transformation matrix satisfies the auxiliary DEs.*)


Print[ And[
  Simplify[D[T,msq]+T . Coefficient[MatrixA[msq],eps,0]] === ConstantArray[0,{nMIs,nMIs}],
  Simplify[D[T,s]+T . Coefficient[MatrixA[s],eps,0]] === ConstantArray[0,{nMIs,nMIs}]
]];


(* ::Text:: *)
(*We define the domain of the variables in order to perform simplifications.*)


Domain = {s<0,msq>0};


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


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


(* ::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]msq = B[msq].g .*)


MatrixB[s] = Table[Coefficient[fder[s][[i]],g],{i,nMIs}];
MatrixB[msq] = Table[Coefficient[fder[msq][[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*A[s].f ,*)
(**)
(*where*)
(**)
(*A[s] = B[s] . Inverse[T] / eps ,*)
(**)
(*and similarly for msq.*)


MatrixA[s] = Simplify[MatrixB[s] . Inverse[T] /eps,Domain];
MatrixA[msq] = Simplify[MatrixB[msq] . Inverse[T]/eps,Domain];

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


(* ::Text:: *)
(*Check the factorisation of eps, namely that MatrixA[s] and MatrixA[msq] are free of eps:*)


Print[ FreeQ[MatrixA[s],eps] && FreeQ[MatrixA[msq],eps] ];


(* ::Text:: *)
(*Alternatively, we may have obtained the new DEs through the "gauge" transformation*)
(**)
(*eps*A'[s] = T . A[s] . T^(-1) + \[PartialD]T/\[PartialD]s . T^(-1)  ,*)
(**)
(*and similarly for msq.*)


(* ::Text:: *)
(*We now check integrability conditions and scaling relation. Note that the integrability conditions*)
(**)
(*eps^2 (A[msq].A[s]-A[s].A[msq]) + eps (D[A[msq],s] - D[A[s],msq]) = 0*)
(**)
(* hold for any value of eps. Thanks to the factorisation of eps, this implies that *)
(**)
(*A[msq].A[s] - A[s].A[msq] = 0*)
(**)
(*and *)
(**)
(*D[A[msq],s] - D[A[s],msq] = 0*)
(**)
(*have to be satisfied separately.*)


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

Scaling = Simplify[s*MatrixA[s]+msq*MatrixA[msq]];
Print["scaling: ", MatrixForm[Scaling]];

Clear[Scaling,IntegrabilityConditions];


(* ::Text:: *)
(*We rewrite the DEs in differential form*)
(**)
(*d f = eps (d Atilde) . f .*)
(**)
(* We check that*)
(**)
(*\[PartialD]Atilde/\[PartialD]s = A[s]*)
(*\[PartialD]Atilde/\[PartialD]msq = A[msq]*)


Atilde = {{- Log[msq], 0},{-2*Log[(Sqrt[1-4*msq/s]-1)/(Sqrt[1-4*msq/s]+1)], - Log[4*msq-s]}};

Print["Atilde = ", MatrixForm[Atilde]];

Print["check Atilde matrix: ",
  And[ 
    Simplify[MatrixA[s]-D[Atilde,s],Domain]===ConstantArray[0,{nMIs,nMIs}], 
    Simplify[MatrixA[msq]-D[Atilde,msq],Domain]===ConstantArray[0,{nMIs,nMIs}]
  ]
];


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


(* ::Text:: *)
(*We change integral basis one more time, from g to f = T . g, with*)


T = msq^eps/Gamma[eps]*{{1,0},{0,Sqrt[(-s)*(4*msq-s)]}};
f = T . g;


(* ::Text:: *)
(*There are two differences with respect to the basis considered in section c. *)
(*First, the factor of 1/Gamma[eps]. Since it is an overall constant prefactor, it is irrelevant for the DEs, and plays a role only for the boundary values. *)
(*Second, the factor of msq^eps. The latter ensures that the new basis is dimensionless. By changing variables from (s,msq) to (msq,x) with*)


Domain = Append[Domain,0<x<1];
s2x = s -> -msq*(1-x)^2/x;


(* ::Text:: *)
(*the basis depends only on the dimensionless variable x. We can thus obtain the connection matrix "Atilde" for this basis by simply changing variables and setting msq = 1 in the connection matrix obtained in section c.*)


Atilde = Simplify[Atilde /. s2x /. msq->1,Domain]/.Log[a_]:>Log[Factor[a]];
Print["Atilde = ", Atilde];


(* ::Chapter:: *)
(*e) Solution of the DE*)


(* ::Text:: *)
(*We spell out the derivative of f w.r.t. x as*)
(**)
(*\[PartialD]f/\[PartialD]x = eps * A'[x] . f ,*)
(**)
(*with *)
(**)
(*A'[x] = D[Atilde,x] .*)
(**)
(*We further perform a partial fraction decomposition to make the pole structure more transparent.*)


MatrixAprime[x] = Apart[D[Atilde,x],x];
Print["A'[x] = ", MatrixAprime[x]];


(* ::Text:: *)
(*We choose the maximum order in eps for the expansion of the solution:*)


MaxEpsOrder = 2;


(* ::Text:: *)
(*We prepare the boundary values at x=x0=1 (s=0) order by order in eps:*)


x0 = 1;
BdryValues[0] = {1,0};

Do[
  BdryValues[w] = {0,0};,
{w,1,MaxEpsOrder}];


(* ::Text:: *)
(*We integrate the DE iteratively order by order in eps, up to the chosen maximum order. *)
(*Assuming that x>1 makes the expressions simpler. The result may then be analytically continued to 0<x<1, if necessary.*)


Clear[f$];

f$[0] = BdryValues[0];

For[w=1, w<=MaxEpsOrder, w++,
  Print["order eps^", w];
  
  f$[w] = Integrate[MatrixAprime[x] . f$[w-1]/.x->t,{t,x0,x},Assumptions->x>1]+BdryValues[w];
 
  Print["solution: ", f$[w]]; 
  Print["------------------------------------------"];
];
  
fAnalytic = Sum[f$[w]*eps^w,{w,0,MaxEpsOrder}];

Clear[f$];


(* ::Text:: *)
(*Finally, our result for the massive bubble integral F_2 is *)


F2=(Inverse[T] . fAnalytic)[[2]];
Print["F_2 = ", F2, "+O[eps^", MaxEpsOrder, "]"];
