(* ::Package:: *)

(* Leading-order psi psi -> psi psi scattering in
     L_int = lambda [(\bar psi psi)^2 - (\bar psi gamma5 psi)^2].

   Requires Mathematica 13.x and FeynCalc 10.x.
   FeynArts is loaded to match the assignment environment.  For this
   four-fermion contact operator, the two tree contributions are the two
   inequivalent Wick contractions of the identical external fermions,
   conventionally labelled t and u.
*)

If[!MemberQ[$Packages, "FeynCalc`"],
   $LoadFeynArts = True;
   Get["FeynCalc`"];
];

$FCDefaultDiracGammaScheme = "NDR";
FCClearScalarProducts[];

ClearAll[
  lambda, p1, p2, p3, p4, s, t, u,
  St, Pt, Su, Pu, Mt, Mu, Mtree,
  trSSt, trPPt, trSPt, trPSt,
  trSSu, trPPu, trSPu, trPSu,
  trSSInt, trSPInt, trPSInt, trPPInt,
  tSpinSumBare, uSpinSumBare, interferenceBare,
  spinSum, spinAverage, finalResult
];

Get[FileNameJoin[{DirectoryName[$InputFileName], "mandelstam_defs.wl"}]];

(* ---------------------------------------------------------------------- *)
(* 1. Tree-level contact diagrams / Wick contractions                     *)
(* ---------------------------------------------------------------------- *)

contactDiagrams = {
  HoldForm[
    "t-contact: (bar u3 Gamma u1)(bar u4 Gamma u2)"
  ],
  HoldForm[
    "u-contact: -(bar u4 Gamma u1)(bar u3 Gamma u2)"
  ]
};

Print["Tree-level contact contractions:"];
Print /@ contactDiagrams;

(* Spinor-level structures.  They are kept as documentation of the
   amplitude; the spin sums below are evaluated directly as Dirac traces. *)
St = HoldForm[(SpinorUBar[p3, 0].SpinorU[p1, 0])*
              (SpinorUBar[p4, 0].SpinorU[p2, 0])];
Pt = HoldForm[(SpinorUBar[p3, 0].GA[5].SpinorU[p1, 0])*
              (SpinorUBar[p4, 0].GA[5].SpinorU[p2, 0])];
Su = HoldForm[(SpinorUBar[p4, 0].SpinorU[p1, 0])*
              (SpinorUBar[p3, 0].SpinorU[p2, 0])];
Pu = HoldForm[(SpinorUBar[p4, 0].GA[5].SpinorU[p1, 0])*
              (SpinorUBar[p3, 0].GA[5].SpinorU[p2, 0])];

Mt = 2 lambda HoldForm[St - Pt];
Mu = -2 lambda HoldForm[Su - Pu];
Mtree = HoldForm[Mt + Mu];

Print["Amplitude structures:"];
Print["Mt = ", Mt];
Print["Mu = ", Mu];
Print["M  = ", Mtree];

(* Helper: evaluate a trace and convert scalar products to s,t,u. *)
toMandelstam[expr_] := Expand[
  (DiracSimplify[DiracTrace[expr], DiracTraceEvaluate -> True] //
     Contract // ExpandScalarProduct) /. mandelstamRules
];

(* ---------------------------------------------------------------------- *)
(* 2. Explicit t-contraction traces                                       *)
(* ---------------------------------------------------------------------- *)

trSSt = toMandelstam[GS[p3].GS[p1]] *
        toMandelstam[GS[p4].GS[p2]];

trPPt = toMandelstam[GS[p3].GA[5].GS[p1].GA[5]] *
        toMandelstam[GS[p4].GA[5].GS[p2].GA[5]];

trSPt = toMandelstam[GS[p3].GS[p1].GA[5]] *
        toMandelstam[GS[p4].GS[p2].GA[5]];

trPSt = toMandelstam[GS[p3].GA[5].GS[p1]] *
        toMandelstam[GS[p4].GA[5].GS[p2]];

Print["t traces:"];
Print["Tr(SS)_t = ", trSSt];
Print["Tr(PP)_t = ", trPPt];
Print["Tr(SP)_t = ", trSPt];
Print["Tr(PS)_t = ", trPSt];

tSpinSumBare = Expand[trSSt + trPPt - trSPt - trPSt];
Print["Sum_spins |St-Pt|^2 = ", tSpinSumBare];

(* ---------------------------------------------------------------------- *)
(* 3. Explicit u-contraction traces                                       *)
(* ---------------------------------------------------------------------- *)

trSSu = toMandelstam[GS[p4].GS[p1]] *
        toMandelstam[GS[p3].GS[p2]];

trPPu = toMandelstam[GS[p4].GA[5].GS[p1].GA[5]] *
        toMandelstam[GS[p3].GA[5].GS[p2].GA[5]];

trSPu = toMandelstam[GS[p4].GS[p1].GA[5]] *
        toMandelstam[GS[p3].GS[p2].GA[5]];

trPSu = toMandelstam[GS[p4].GA[5].GS[p1]] *
        toMandelstam[GS[p3].GA[5].GS[p2]];

Print["u traces:"];
Print["Tr(SS)_u = ", trSSu];
Print["Tr(PP)_u = ", trPPu];
Print["Tr(SP)_u = ", trSPu];
Print["Tr(PS)_u = ", trPSu];

uSpinSumBare = Expand[trSSu + trPPu - trSPu - trPSu];
Print["Sum_spins |Su-Pu|^2 = ", uSpinSumBare];

(* ---------------------------------------------------------------------- *)
(* 4. Explicit t-u interference traces                                    *)
(* ---------------------------------------------------------------------- *)

trSSInt = toMandelstam[
  GS[p3].GS[p1].GS[p4].GS[p2]
];

trSPInt = toMandelstam[
  GS[p3].GS[p1].GA[5].GS[p4].GS[p2].GA[5]
];

trPSInt = toMandelstam[
  GS[p3].GA[5].GS[p1].GS[p4].GA[5].GS[p2]
];

trPPInt = toMandelstam[
  GS[p3].GA[5].GS[p1].GA[5].GS[p4].GA[5].GS[p2].GA[5]
];

Print["t-u interference traces:"];
Print["T_SS = ", trSSInt];
Print["T_SP = ", trSPInt];
Print["T_PS = ", trPSInt];
Print["T_PP = ", trPPInt];

(* The scalar/pseudoscalar combination gives T_SS-T_SP-T_PS+T_PP = 0. *)
interferenceBare = Expand[
  trSSInt - trSPInt - trPSInt + trPPInt
];
Print["Interference trace combination = ", interferenceBare];

(* ---------------------------------------------------------------------- *)
(* 5. Couplings, identical-fermion sign, and initial-spin average         *)
(* ---------------------------------------------------------------------- *)

(* M = 2 lambda [(St-Pt) - (Su-Pu)].
   The interference term is zero, but its coefficient is retained explicitly. *)
spinSum = Expand[
  (2 lambda)^2 tSpinSumBare +
  (2 lambda)^2 uSpinSumBare -
  2 (2 lambda)^2 interferenceBare
];

spinAverage = Expand[spinSum/4];

(* Minimal massless polynomial. *)
finalResult = Factor[
  spinAverage /. u^2 -> u^2
];

Print["Spin-summed |M|^2 = ", spinSum];
Print["Spin-averaged |M|^2 = ", spinAverage];
Print["Final result = ", finalResult];

expectedResult = 8 lambda^2 (t^2 + u^2);
If[Expand[finalResult - expectedResult] =!= 0,
  Print["ERROR: result does not match the expected polynomial."];
  Abort[];
];

Put[
  finalResult,
  FileNameJoin[{DirectoryName[$InputFileName], "amplitude_squared.wl"}]
];

finalResult