QED vacuum polarization to two loops

Needs["HighEnergyPhysics`FeynCalc`"];

FeynCalc 4.1.1    For help, type ?FeynCalc,  use the built-in help system  or visit www.feyncalc.org

The expression we wish to consider is the product of the two factors below:

tra = 
      DiracTrace[
          GAD[ν,  σ,  μ,  λ,  α,  ρ,  β,  λ]] // FCI

notrace = -e^4 FAD[q, p, p + k, p + q + k, p + q] FVD[p, β] 
            FVD[p + k, α] FVD[p + q + k, μ] 
            FVD[p + q, ν] /. {α -> ν, β -> μ, ν -> α, μ -> β} (*Seems to be a bug in Itzykson & Zuber*)  // FCI

The trace is easily evaluated:

trace = 
      Collect[Simplify [
                Tr[tra]] //. {(6  a_ - D a_) -> (D - 6) 
                    min a, (D a_ - 6  a_) -> (D - 6) 
                    a, (2  a_ - D a_) -> (D - 2) 
                    min a, (D a_ - 2  a_) -> (D - 2) 
                    a}, {(D - 2), (D - 6)}] /. min -> -1 // 
        Simplify

The other factor will be Feynman parameterized. First, the structure of the expression:

(FeynmanParametrize1[notrace, p, Method -> Exp, 
                Integrate -> False, CompleteSquare -> False, 
                FeynmanParameterNames -> {α1, α2, α3, α4}] // ExpandScalarProduct // 
            Simplify) /. {α1 + α2 + α3 + α4 -> Σ, -α1 - α2 - α3 - α4 -> -Σ} // Simplify

(FeynmanParametrize1[notrace, p, Method -> Exp, 
                Integrate -> False, CompleteSquare -> True, 
                FeynmanParameterNames -> {α1, α2, α3, α4}] // ExpandScalarProduct // 
            Simplify) /. {α1 + α2 + α3 + α4 -> Σ, -α1 - α2 - α3 - α4 -> -Σ} // Simplify

Tell FeynmanParametrize1 to do the p integral:

notrace1 = 
(FeynmanParametrize1[notrace, p, 
                Method -> Exp, Integrate -> True, 
                FeynmanParameterNames -> {α1, α2, α3, α4}] // ExpandScalarProduct // Simplify);

Multiply on the trace and contract Lorentz indices:

qexpr = 
        notrace1 /. (a : Dot[Integratedx[__], __]) :> 
            ReplacePart[a, 
              Contract[a[[-1]] trace], -1];

Simplify a bit:

qexprcc = 
        qexpr /. (aa : Dot[Integratedx[__], __]) :> 
            MapAt[((Collect[# /.  
                            2^(a_ - D) -> (2^a) tD, {tD, 
                            e, E^_, π^_, _^(-(D/2) + _), Pair[__]}] //. p : (Pair[__]*_) :> Simplify[p] /. 
                      tD -> 2^(-D)) &), 
              Expand[ExpandScalarProduct[aa]], -1];

qexprcc

Tell FeynmanParametrize1 to do the q integral:

pexpr = 
        FeynmanParametrize1[qexprcc, q, Method -> Exp, 
          Flatten -> True, Integrate -> True, 
          FeynmanParameterNames -> {x}];

Simplify a bit:

pexpr // LeafCount

tmp = pexpr //. {p : (Pair[a_LorentzIndex, b_] 
                      Pair[a_LorentzIndex, c_]) :> Contract[p], 
              p : (Pair[_LorentzIndex, _]^2) :> 
                Contract[p]} /. 
          E^a_ :> E^Simplify[a];

tmp // LeafCount

pexprcc = 
        tmp /. (aa : 
                Dot[Integratedx[__], ___, _?(FreeQ[#, 
                          Integratedx] &)]) :> 
            MapAt[((Collect[# /.  
                              2^(a_ - 2  D) :> (2^
                                      a) t2D, {E^_, e, 
                              t2D, π^_, _Pair, _^(_?(! FreeQ[#, 
                                        D] &))}] /. 
                          p : (_Pair*_) :> Simplify[p] /. 
                        t2D -> 2^(-2 D) /. 
                      a : ((_?(! FreeQ[#, 
                                        x] &))^_) :> 
                        Simplify[a]) &), aa, -1];

pexprcc // LeafCount

pexps = pexprcc // Simplify;

pexpss = 
        pexps /. p : (a_Pair*b_) :> 
            a*FullSimplify[b];

pexpss // LeafCount

Use FeynmanReduce to get rid of the exponentials and shrink the integration region:

rexp = 
      FeynmanReduce[
        pexpss, {α1, α2, α3, α4}]

Variable substitutions and integration: ***Takes very long***

\(\(fdi = 
        CheckF[FeynmanIntegrate[
            rexp /. {α3 -> α4, α4 -> α3}, {k}, {α1, α2, α3, α4, 
              x}], "TwoLoopVacPolfdi"];\)\)

sdi = Series[fdi /. _Dot -> 0, {Epsilon, 0, 0}] // Normal // 
        FullSimplify