(* ::Package:: *)

(* ::Subtitle:: *)
(*Restriction z1 \[RightArrow] 0 in the massless limit*)


SetDirectory @ NotebookDirectory[];
Get["../utils.m"]


(* ::Section:: *)
(*Pfaffian matrices for the massive Bhabha double box*)


{P1, P2} = Import["massive_pfaffians.m"];


(* ::Text:: *)
(*P1 not in normal form: 3rd order pole (want 1st order)*)


P1 // mseries1[z1, 0]


(* ::Text:: *)
(*P2 not in normal form: 2nd order pole (want finite)*)


P2 // mseries1[z1, 0]


(* ::Section:: *)
(*Normal form by Moser reduction*)


(* ::Subsection:: *)
(*Shearing gauge transformations*)


ClearAll[G1, G2, G3];
G1 = DiagonalMatrix[Join[
    ConstantArray[z1, 8],
    {}
] // PadRight[#, 23, 1]&];
G2 = DiagonalMatrix[Join[
    ConstantArray[z1, 12],
    {}
] // PadRight[#, 23, 1]&];
G3[1] = DiagonalMatrix[Join[
    ConstantArray[z1^-1, 3],
    {}
] // PadRight[#, 23, 1]&];
G3[2] = DiagonalMatrix[Join[
    ConstantArray[z1^-1, 5],
    {}
] // PadRight[#, 23, 1]&];
G3[3] = DiagonalMatrix[Join[
    ConstantArray[1, 8],
    ConstantArray[z1^-1, 2],
    {}
] // PadRight[#, 23, 1]&];
G3[4] = DiagonalMatrix[Join[
    ConstantArray[1, 9],
    ConstantArray[z1^-1, 4],
    {}
] // PadRight[#, 23, 1]&];
G3[5] = DiagonalMatrix[Join[
    ConstantArray[1, 8],
    ConstantArray[z1^-1, 1],
    {}
] // PadRight[#, 23, 1]&];
G3[6] = DiagonalMatrix[Join[
    ConstantArray[1, 8],
    ConstantArray[z1^-1, 3],
    {}
] // PadRight[#, 23, 1]&];
G3[7] = DiagonalMatrix[Join[
    ConstantArray[1, 8],
    ConstantArray[z1^-1, 6],
    {}
] // PadRight[#, 23, 1]&];


(* ::Subsection:: *)
(*Remove higher order poles*)


(* ::Text:: *)
(*Remove 3rd order pole in P1*)


P1G1 = P1 // gauge[G1, z1] // Factor;


P1G1 // mseries1[z1, 0]


(* ::Text:: *)
(*Remove 2nd order pole in P1*)


P1G2 = P1G1 // gauge[G2, z1] // Factor;


P1G2 // mseries1[z1, 0]


(* ::Subsection:: *)
(*Remove the resonances in the spectrum of P1*)


P1G2 // res[z1] // JordanDecomposition // Last // Diagonal


ClearAll[P1step, P1jord, JordanG];


(* ::Text:: *)
(*Look at the eigenvalues of P1 and derive the required gauge transformation G3*)


JordanG[1] = P1G2 // res[z1] // JordanDecomposition // First;
P1jord[1] = P1G2 // conj[JordanG[1]] // Factor;
P1step[1] = P1jord[1] // gauge[G3[1], z1] // Factor;


(* ::Text:: *)
(*We just reduced the -2 eigenvalues down to -1*)


P1step[1] // res[z1] // JordanDecomposition // Last // Diagonal


(* ::Subsubsection:: *)
(*Repeat*)


JordanG[2] = P1step[1] // res[z1] // JordanDecomposition // First;
P1jord[2] = P1step[1]  // conj[JordanG[2]] // Factor;
P1step[2] = P1jord[2] // gauge[G3[2], z1] // Factor;


JordanG[3] = P1step[2] // res[z1] // JordanDecomposition // First;
P1jord[3] = P1step[2]  // conj[JordanG[3]] // Factor;
P1step[3] = P1jord[3] // gauge[G3[3], z1] // Factor;


JordanG[4] = P1step[3] // res[z1] // JordanDecomposition // First;
P1jord[4] = P1step[3]  // conj[JordanG[4]] // Factor;
P1step[4] = P1jord[4] // gauge[G3[4], z1] // Factor;


JordanG[5] = P1step[4] // res[z1] // JordanDecomposition // First;
P1jord[5] = P1step[4]  // conj[JordanG[5]] // Factor;
P1step[5] = P1jord[5] // gauge[G3[5], z1] // Factor;


JordanG[6] = P1step[5] // res[z1] // JordanDecomposition // First;
P1jord[6] = P1step[5]  // conj[JordanG[6]] // Factor;
P1step[6] = P1jord[6] // gauge[G3[6], z1] // Factor;


JordanG[7] = P1step[6] // res[z1] // JordanDecomposition // First;
P1jord[7] = P1step[6]  // conj[JordanG[7]] // Factor;
P1step[7] = P1jord[7] // gauge[G3[7], z1] // Factor;


(* ::Subsubsection:: *)
(*Result*)


P1G = P1step[7];


(* ::Text:: *)
(*Still have the 1st order pole*)


P1G // mseries1[z1, 0] 


(* ::Text:: *)
(*The spectrum is now free of resonances*)


P1G // res[z1] // Eigenvalues


(* ::Text:: *)
(*The list of required gauge transformations*)


Gs = List[
    G1,
    G2,
    JordanG[1], G3[1],
    JordanG[2], G3[2],
    JordanG[3], G3[3],
    JordanG[4], G3[4],
    JordanG[5], G3[5],
    JordanG[6], G3[6],
    JordanG[7], G3[7],
    Nothing
] // Map[SparseArray];


(* ::Subsection:: *)
(*Apply gauge transformations to P2 and save*)


P2G = Fold[gauge[#2, z2][#1]& /* Normal /* Factor /* TM["gauge"], P2, Gs];


Export["pfaffian_physical_2L_bhabha_normal.m", {P1G, P2G}]
Export["pfaffian_physical_2L_bhabha_gauge_transformations.m", Gs]


(* ::Section:: *)
(*Restriction*)


(* ::Text:: *)
(*Import saved results from above*)


{P1G, P2G} = Import["massive_pfaffians_normal_form.m"];


Gs = Import["gauge_transformations.m"];


(* ::Text:: *)
(*Combine the gauge transfromations*)


ClearAll @ G;
Ginv = Fold[Inverse[#2] . #1&, IdentityMatrix[23], Gs];
G = Fold[Dot, IdentityMatrix[23], Gs];


(* ::Text:: *)
(*Master integrals for the massive system*)


mastersm = Import["massive_master_integrals.m"]


(* ::Text:: *)
(*Master integrals for the massless system*)


masters0 = Import["massless_master_integrals.m"]


(* ::Text:: *)
(*masters0 = B . mastersm*)


B = CoefficientArrays[masters0, mastersm] [[2]] // Normal;
B // MatrixForm


BG = B . G // mseries1[z1, 0] // Part[#, 2, 1]&;


(* ::Text:: *)
(*R-matrix*)


R = P1G // res[z1] // rowReduce


(* ::Text:: *)
(*M-matrix*)


M = Join[BG, R]


(* ::Text:: *)
(*Restriction*)


Q2 = (D[M,z2] + M . (P2G /. z1 -> 0)) . Inverse[M] // RightComposition[
    #[[;;8 , ;;8]]&,
    Factor,
    Identity
];


Q2answer = Import["massless_pfaffians_answer.m"];


Q2 - Q2answer // Factor
