(* ::Package:: *)

Quit


(* s=(p1+p2)^2, t=(p1-p3)^2, u=(p2-p3)^2=-s-t *)


(* ::Subsection::Closed:: *)
(*Setup reduction*)


SetDirectory["~/applications/fire/FIRE5"];
Get["FIRE5.m"];
Internal={l};
External={p1,p2,p3};
Propagators={l^2,(l+p1)^2,(l+p1+p2)^2,(l+p3)^2};
Replacements={p1^2->0,p2^2->0,p3^2->0,p1*p2->s/2,p1*p3->-t/2,p2*p3->(s+t)/2};
PrepareIBP[];
Prepare[AutoDetectRestrictions->True];
SaveStart["~/Projects/MultiLoop/Examples/MasslessBox"];
Quit[];


(* ::Subsection::Closed:: *)
(*Differential equation and canonical basis*)


SetDirectory["~/applications/fire/FIRE5"];
Get["FIRE5.m"];
LoadStart["~/Projects/MultiLoop/Examples/MasslessBox",1];
Burn[];


F[1,{2,2,2,2}]


MasterIntegrals[]


(* Choose dimensionless basis integrals by rescaling with (-s)^dim s.t. they only depend on the dimensionless ratio x=t/s *)
(* f1 = (-s)^eps*F[1,{0,1,0,1}], f2 = (-s)^eps*F[1,{1,0,1,0}], f3 = (-s)^(2+eps)*F[1,{1,1,1,1}] *)


(* Build derivative operator Dx *)

(* scalar product *)
Attributes[S] = {Orderless};
S /: S[a_,{mu}]^2 := S[a,a];
S /: S[a_,{mu}]*S[b_,{mu}] := S[a,b];
S[a_Plus,b_] := Sum[S[a[[i]],b],{i,1,Length[a]}];
S[n_?NumericQ*a_,b_] := n*S[a,b];

(* derivatives of scalar producs *)
diff = {\!\(\*SuperscriptBox[\(S\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[p1,k_]->S[{mu},k],\!\(\*SuperscriptBox[\(S\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[k_,p1]->S[k,{mu}]};

(* express scalar products through kinematic invariants or denominators *)
reps = {S[p1,p1]->0,S[p2,p2]->0,S[p3,p3]->0,S[p1,p2]->s/2,S[p1,p3]->-t/2,S[p2,p3]->(s+t)/2};
todens = {S[l,l]->D1,S[l,p1]->(D2-D1)/2,S[l,p2]->(D3-D2-s)/2,S[l,p3]->(D4-D1)/2};

(* Determine coefficients in ansatz Dx[...] = (a*S[p1,{mu}]+b*S[p2,{mu}]+c*S[p3,{mu}])*D[...,p1] *)
Solve[
(* Dx[x] \[Equal] 1 *)
1 == FullSimplify[Expand[(a*S[p1,{mu}]+b*S[p2,{mu}]+c*S[p3,{mu}])*D[S[p1-p3,p1-p3]/(S[p1,p1]+2S[p1,p2]+S[p2,p2]),p1] /. diff] /. reps /. {t->x*s}] &&
(* Dx[p1^2] \[Equal] 0 *)
0 == FullSimplify[Expand[(a*S[p1,{mu}]+b*S[p2,{mu}]+c*S[p3,{mu}])*D[S[p1,p1],p1] /. diff] /. reps /. {t->x*s}] &&
(* Dx[p4^2] \[Equal] 0 *)
0 == FullSimplify[Expand[(a*S[p1,{mu}]+b*S[p2,{mu}]+c*S[p3,{mu}])*D[S[p1+p2-p3,p1+p2-p3],p1] /. diff] /. reps /. {t->x*s}],
{a,b,c}]


(* The derivative wrt x takes the form *)
Dx[expr_] := FullSimplify[Expand[((1-x)/(2*x*(1+x))*S[p1,{mu}] - 1/(2*(1+x))*S[p2,{mu}] - 1/(2*x*(1+x))*S[p3,{mu}])*D[expr,p1] /. diff] /. reps /. {t->x*s}]


(* Project FIRE's master integrals on our basis choice *)
ToBasisF = {G[1,{0,1,0,1}]->(-s)^-eps*f1,
   G[1,{1,0,1,0}]->(-s)^-eps*f2,
   G[1,{1,1,1,1}]->(-s)^(-2-eps)*f3};


(* d/(dx) {f1,f2,f3} *)
A = Apart[Expand[Dx[{
   (-S[p1+p2,p1+p2])^eps/(S[l+p1,l+p1]*S[l+p3,l+p3]),
   (-S[p1+p2,p1+p2])^eps/(S[l,l]*S[l+p2,l+p2]),
   (-S[p1+p2,p1+p2])^(2+eps)/(S[l,l]*S[l+p1,l+p1]*S[l+p1+p2,l+p1+p2]*S[l+p3,l+p3])}
   ]*Faux[1,{0,0,0,0}] /. todens] //. {
      Faux[1,{a_,b_,c_,d_}]*D1 -> Faux[1,{a-1,b,c,d}],Faux[1,{a_,b_,c_,d_}]*D1^n_ -> Faux[1,{a-n,b,c,d}],
      Faux[1,{a_,b_,c_,d_}]*D2 -> Faux[1,{a,b-1,c,d}],Faux[1,{a_,b_,c_,d_}]*D2^n_ -> Faux[1,{a,b-n,c,d}],
      Faux[1,{a_,b_,c_,d_}]*D3 -> Faux[1,{a,b,c-1,d}],Faux[1,{a_,b_,c_,d_}]*D3^n_ -> Faux[1,{a,b,c-n,d}],
      Faux[1,{a_,b_,c_,d_}]*D4 -> Faux[1,{a,b,c,d-1}],Faux[1,{a_,b_,c_,d_}]*D4^n_ -> Faux[1,{a,b,c,d-n}]
   } /. {Faux->F} /. {d->4-2*eps,t->x*s} /. ToBasisF /. {
   {a_,b_,c_} :> {Coefficient[a,{f1,f2,f3}],Coefficient[b,{f1,f2,f3}],Coefficient[c,{f1,f2,f3}]}},
x];
A//MatrixForm


(* The first step is to bring the basis into dlog form, to do so we only have to remove the 1/x^2 term *)
(* For a basis transformation fprime = Tg, the matrix A changes as follows *)
Aprime[A_,T_] := Apart[Inverse[T].A.T - Inverse[T].D[T,x],x];


(* This implies the condition -T.Aprime+A.T-Dx[T] \[Equal] 0 where Aprime = a/x + b/(1+x) *)
(* Since we want to remove the term with 1/x^2 we try the simple ansatz diag(x^t1,x^t2,x^t3) *)
With[{T = {{x^t1,0,0},{0,x^t2,0},{0,0,x^t3}}, 
Aprime = Table[ToExpression["a"<>ToString[i]<>ToString[j]]/x + ToExpression["b"<>ToString[i]<>ToString[j]]/(1+x),{i,1,3},{j,1,3}]},
Collect[-T.Aprime + A.T - D[T,x],x]]


(* This suggests the following form *)
Aprime[A,{{1,0,0},{0,1,0},{0,0,1/x}}]


(* To transform this to canonical form we need a further transformation that only depends on eps. 
To make the off-diagonal entries O(eps) we use *)
Aprime[A,{{1,0,0},{0,1,0},{0,0,1/(x*eps)}}]
(* Then to get rid of the (1-2*eps) factors *)
Aprime[A,{{-1/(1-2*eps),0,0},{0,-1/(1-2*eps),0},{0,0,1/(x*eps)}}]


(* Check new basis. Extra eps*Exp[eps*EulerGamma] to have eps-expansion start at finite order and to remove EulerGamma terms. This does not change DE. *)
ToBasisG = {G[1,{0,1,0,1}]->(g1 (-s)^-eps)/((-1+2 eps)eps*Exp[eps*EulerGamma]),
   G[1,{1,0,1,0}]->(g2 (-s)^-eps)/((-1+2 eps)eps*Exp[eps*EulerGamma]),
   G[1,{1,1,1,1}]->(g3 (-s)^-eps)/(eps^2*Exp[eps*EulerGamma]*s^2*x)};

(* Generate diff eq for new basis *)
Apr = Apart[Expand[Dx[eps*Exp[eps*EulerGamma]*{((-S[p1+p2,p1+p2])^eps*S[p1-p3,p1-p3])/(S[l+p1,l+p1]*S[l+p3,l+p3]^2),(-S[p1+p2,p1+p2])^eps/(S[l,l]*S[l+p2,l+p2]^2),(eps*(-S[p1+p2,p1+p2])^eps*S[p1+p2,p1+p2]*S[p1-p3,p1-p3])/(S[l,l]*S[l+p1,l+p1]*S[l+p1+p2,l+p1+p2]*S[l+p3,l+p3])}]*Faux[1,{0,0,0,0}] /. todens] //. {
   Faux[1,{a_,b_,c_,d_}]*D1 -> Faux[1,{a-1,b,c,d}],Faux[1,{a_,b_,c_,d_}]*D1^n_ -> Faux[1,{a-n,b,c,d}],
   Faux[1,{a_,b_,c_,d_}]*D2 -> Faux[1,{a,b-1,c,d}],Faux[1,{a_,b_,c_,d_}]*D2^n_ -> Faux[1,{a,b-n,c,d}],
   Faux[1,{a_,b_,c_,d_}]*D3 -> Faux[1,{a,b,c-1,d}],Faux[1,{a_,b_,c_,d_}]*D3^n_ -> Faux[1,{a,b,c-n,d}],
   Faux[1,{a_,b_,c_,d_}]*D4 -> Faux[1,{a,b,c,d-1}],Faux[1,{a_,b_,c_,d_}]*D4^n_ -> Faux[1,{a,b,c,d-n}]} /. {Faux->F} /. {d->4-2*eps,t->x*s} /. ToBasisG /. {
   {a_,b_,c_} :> {Coefficient[a,{g1,g2,g3}],Coefficient[b,{g1,g2,g3}],Coefficient[c,{g1,g2,g3}]}},x];
Apr//MatrixForm


(* Agrees with result from transformation *)
Aprime[A,{{-1/(1-2*eps),0,0},{0,-1/(1-2*eps),0},{0,0,1/(x*eps)}}] - Apr


(* Coefficients of eps/x and eps/(1+x) *)
a = Limit[Coefficient[Coefficient[Apr,eps],x,-1],x->Infinity]
b = Coefficient[Coefficient[Apr,eps],1+x,-1]


(* ::Subsection::Closed:: *)
(*Solution*)


(* Differential equation does not provide any information on g1,g2 since they are single-scale integrals. Result from Feynman parameters *)
g1 = -x^(-eps)*Exp[eps*EulerGamma]*Gamma[1-eps]^2*Gamma[1+eps]/Gamma[1-2*eps]; 
g2 = -Exp[eps*EulerGamma]*Gamma[1-eps]^2*Gamma[1+eps]/Gamma[1-2*eps];
g1exp = FullSimplify[Series[g1,{eps,0,3}]]
g2exp = FullSimplify[Series[g2,{eps,0,3}]]


(* From the absence of u-channel cuts we obtain the boundary condition g3[-1] = -2(g1[-1]+g2[-1]) to all orders in eps *)
g3BC = FullSimplify[Series[-2g1-2g2,{eps,0,3}]/.{Log[x]->I \[Pi]}]


(* See notes for solution in terms of Goncharov polylogs *)


(* Express G[-1,0,0;x] in terms of Log and PolyLog *)
Integrate[1/(1+t)*Log[t]^2/2,{t,0,x},Assumptions->x>0]


(* Determine boundary condition at 0 from boundary condition at -1 *)
Solve[Limit[((7 \[Pi]^2)/6 Log[x]+Log[x]^3/3-(\[Pi]^2+Log[x]^2)*Log[1+x]-2 Log[x] PolyLog[2,-x]+2 PolyLog[3,-x]+g330),x->-1] == Coefficient[g3BC,eps,3],g330]


(* Result for g3 up to eps^3 *)
4 - 2*Log[x]*eps - (4 \[Pi]^2)/3*eps^2 + ((7 \[Pi]^2)/6 Log[x]+Log[x]^3/3-(\[Pi]^2+Log[x]^2)*Log[1+x]-2 Log[x] PolyLog[2,-x]+2 PolyLog[3,-x]-(34 Zeta[3])/3) eps^3