Combined virtual contributions from virtualt_clean.nb, virtualu_clean.nb, virtualtu_clean.nb and virtualself_clean.nb

Original mathematica files are total_virtual_rot_nondynt.nb, (also u, tu) and total_virtual_rot_nondyn_self.nb, total_virtual_rot_nondyn_selfa.nb (both give same self energy contribution)

Overall constants and other terms

factorLO  (* 1/(32pi s) P'/P *) = (( Pp)/(128 π (H^2) P)) (* Note Pp depends on mf *)

factorNLOVirtual = (* 1/(32πs)(P'/P) 1/(2(2pi)^3) is the common factor pulled out of all NLO terms, factor of 2 from 1/(|2k0|) *)  1/(128 (H^2))(Pp/P) 1/(16 (π^4))

factorNLO3 = factorNLOVirtual  (* for dsigmaLO/dmf^2 in self energy B3, F3 to match with other NLO terms, defining dsigma/dm2 w.r.t. P' as (-1/(2P'^2) (A + 3 B P'^2) where sigmaLO = factorLO (A + B P'^2) *)

sigmaDiracLO = InttLO  (* X factorLO *)

sigmaMajoranaLO = InttLO + IntuLO - InttuLO   (* X factorLO *)// Simplify

dsigmaDiracm2 =  -(( 1)/(2(Pp^2))) ( Coefficient[sigmaDiracLO, Pp,0] + 3 (Pp^2)  Coefficient[sigmaDiracLO, Pp,2] ) // Simplify (* with factorLO removed and so get same coefficient of factorNLOVirtual *)

dsigmaMajoranam2 =(* with factorLO removed  and so get same coefficient of factorNLOVirtual *)-(( 1)/(2(Pp^2))) ( Coefficient[sigmaMajoranaLO, Pp,0] + 3 (Pp^2)  Coefficient[sigmaMajoranaLO, Pp,2] ) // Simplify

G part: thermal photon

t-channel part

G1t=((256 (e^2) π (3 (H^4)+(P^2) (Pp^2)) (λ^4))/(3 (mh^6)))

G2t=((64 (e^2) (3 (H^2)+(P^2)) π (λ^4))/(3 (mh^4)))

G3t=((128 (e^2) π (3 (H^4)+(H^2) (3 (mc^2)+4 (P^2))+(P^2) (-(mf^2)+(Pp^2))) (λ^4))/(3 (mh^6)))

G4t=(1/(15 (mf^2) (mh^8)))512 (e^2) π (30 (H^8)-(mf^2) (P^4) (Pp^2)+5 (mc^2) (P^2) (Pp^4)-15 (H^6) ((mc^2)-(mf^2)+2 (Pp^2))+15 (H^4) (3 (mf^2) (P^2)+((mc^2)+2 (P^2)) (Pp^2))+(H^2) (P^2) ((mf^2) (4 (P^2)-15 (Pp^2))-5 (Pp^2) ((mc^2)+6 (Pp^2)))) (λ^4)

G5t=-((512 (e^2) π (3 (H^4)+(P^2) (Pp^2)) (λ^4))/(3 (mh^6)))

u-channel part

G1u=((256 (e^2) π (3 (H^4)+(P^2) (Pp^2)) (λ^4))/(3 (mh^6)))

G2u=((64 (e^2) (3 (H^2)+(P^2)) π (λ^4))/(3 (mh^4)))

G3u=((128 (e^2) π (3 (H^4)+(H^2) (3 (mc^2)+4 (P^2))+(P^2) (-(mf^2)+(Pp^2))) (λ^4))/(3 (mh^6)))

G4u=(1/(15 (mf^2) (mh^8)))512 (e^2) π (30 (H^8)-(mf^2) (P^4) (Pp^2)+5 (mc^2) (P^2) (Pp^4)-15 (H^6) ((mc^2)-(mf^2)+2 (Pp^2))+15 (H^4) (3 (mf^2) (P^2)+((mc^2)+2 (P^2)) (Pp^2))+(H^2) (P^2) ((mf^2) (4 (P^2)-15 (Pp^2))-5 (Pp^2) ((mc^2)+6 (Pp^2)))) (λ^4)

G5u=-((512 (e^2) π (3 (H^4)+(P^2) (Pp^2)) (λ^4))/(3 (mh^6)))

tu cross-channel part

G1tu=((128 (e^2) (mc^2) π ((H^2)+(mf^2)+3 (Pp^2)) (λ^4))/(mh^6))

G2tu=((128 (e^2) (mc^2) π (λ^4))/(mh^4))

G3tu=((128 (e^2) (mc^2) π (3 (H^2)-(mf^2)+(Pp^2)) (λ^4))/(mh^6))

G4tu=(1/(3 (mf^2) (mh^8)))512 (e^2) (mc^2) π ((H^2)+(Pp^2)) (6 (H^4)+(mf^2) (P^2)+3 (mc^2) (Pp^2)-3 (H^2) ((mc^2)-(mf^2)+2 (Pp^2))) (λ^4)

G5tu=-((512 (e^2) (mc^2) π ((H^2)+(Pp^2)) (λ^4))/(mh^6))

G part: thermal fermion

t-channel part

G1ft=-((128 (e^2) (-3 (H^2) (mc^2)+(-4 (H^2)+(mf^2)) (P^2)) π (λ^4))/(3 (mh^6)))

G2ft=-((32 (e^2) (2 (H^2)+(ot^2)) (3 (H^2)+(P^2)) π (λ^4))/(3 (mh^4) (H-ot) (H+ot)))

G2at=-((32 (e^2) (2 (H^2)+(ot^2)) (3 (H^2)+(P^2)) π (λ^4))/(3 (mh^4) (H-ot) (H+ot)))

G3at=-((128 (e^2) (-3 (H^2) (mc^2)+(-4 (H^2)+(mf^2)) (P^2)) π (λ^4))/(3 (mh^6)))

u-channel part

G1fu=((128 (e^2) (-((mf^2) (P^2))+(H^2) (3 (mc^2)+4 (P^2))) π (λ^4))/(3 (mh^6)))

G2fu=-((32 (e^2) (2 (H^2)+(ot^2)) (3 (H^2)+(P^2)) π (λ^4))/(3 (mh^4) (H-ot) (H+ot)))

G2au=-((32 (e^2) (2 (H^2)+(ot^2)) (3 (H^2)+(P^2)) π (λ^4))/(3 (mh^4) (H-ot) (H+ot)))

G3au=((128 (e^2) (-((mf^2) (P^2))+(H^2) (3 (mc^2)+4 (P^2))) π (λ^4))/(3 (mh^6)))

tu cross-channel part

G1ftu=((128 (e^2) (mc^2) (2 (H^2)-(mf^2)) π (λ^4))/(mh^6))

G2ftu=-((64 (e^2) (mc^2) (2 (H^2)+(ot^2)) π (λ^4))/((mh^4) (H-ot) (H+ot)))

G2atu=-((64 (e^2) (mc^2) (2 (H^2)+(ot^2)) π (λ^4))/((mh^4) (H-ot) (H+ot)))

G3atu=((128 (e^2) (mc^2) (2 (H^2)-(mf^2)) π (λ^4))/(mh^6))

K part: needed for soft IR divergence, thermal photon only. These appear to be vanishing but only because the soft and collinear divergences mix and we are considering only the soft here.

t-channel part

K1t=0

K2t=0

K3t=0

K4t=-(1/(15 (mh^8)))512 (e^2) (15 (H^6)-5 (mf^4) (P^2)+(mf^2) (P^4)+15 (H^4) ((mf^2)+2 (P^2))+(H^2) (20 (mf^2) (P^2)+3 (P^4))) π (λ^4)

K5t=0

u-channel part

K1u=ConditionalExpression[0,(((Im[H])/(Im[Pp]))≥1||Re[H]≥((Im[H] Re[Pp])/(Im[Pp]))||((Im[H])/(Im[Pp]))≤-1)&&((H/Pp)∉(Reals)||Re[(H/Pp)]<-1||Re[(H/Pp)]>1)]

K2u=0

K3u=0

K4u=-(1/(15 (mh^8)))512 (e^2) (15 (H^6)-5 (mf^4) (P^2)+(mf^2) (P^4)+15 (H^4) ((mf^2)+2 (P^2))+(H^2) (20 (mf^2) (P^2)+3 (P^4))) π (λ^4)

K5u=0

tu cross-channel part

K1tu=ConditionalExpression[0,(((Im[H])/(Im[Pp]))≥1||Re[H]≥((Im[H] Re[Pp])/(Im[Pp]))||((Im[H])/(Im[Pp]))≤-1)&&((H/Pp)∉(Reals)||Re[(H/Pp)]<-1||Re[(H/Pp)]>1)]

K2tu=0

K3tu=0

K4tu=((512 (e^2) (mc^2) (-2 (H^2)+(mf^2)) (3 (H^2)+3 (mf^2)+(P^2)) π (λ^4))/(3 (mh^8)))

K5tu=0

Self energy part: contributions to both G, K, and Gtot (coll div) from the self-energy parts:

Self Energy Part for 11, 44, 14 (t, u, tu) terms of real photon case. Since same for antifermion self-energy, also for 33, 66, and 16 contributions

Self energy part, finite; already integrated over ω or ot (from both B1, B3, and F1, F3; see virtualself_clean.nb)

totBfinite=ConditionalExpression[((16 (e^2) (π^3) (-6 (H^4)+3 (H^2) (mc^2)+(3 (mc^2)-2 (P^2)) (Pp^2)) (λ^4))/(9 (beta^2) (mh^4) (Pp^2))),Re[beta]>0]

totFfinite=ConditionalExpression[((8 (e^2) (π^3) (-6 (H^4)+3 (H^2) (mc^2)+(3 (mc^2)-2 (P^2)) (Pp^2)) (λ^4))/(9 (beta^2) (mh^4) (Pp^2))),Re[beta]>0]

totBDirac=ConditionalExpression[-((16 (e^2) (π^3) (3 (H^4)+(P^2) (Pp^2)) (λ^4))/(9 (beta^2) (mh^4) (Pp^2))),Re[beta]>0]

totFDirac=ConditionalExpression[-((8 (e^2) (π^3) (3 (H^4)+(P^2) (Pp^2)) (λ^4))/(9 (beta^2) (mh^4) (Pp^2))),Re[beta]>0]

Self energy part, Divergent (from B2 and F2 only)

KBt=((64 (e^2) (mf^2) π (3 (H^4)+(P^2) (Pp^2)) (λ^4))/(3 (mh^4) (H-Pp) (H+Pp) (ω^2)))

(((2 (H^2) - (mf^2))/(4(H^2)) (KBt/. {Pp->(H - (mf^2)/(2H)) })) + O[mf]^4) /. {H-> mc, P-> 0}// Simplify

TeXForm[%]

KBu=((64 (e^2) (mf^2) π (3 (H^4)+(P^2) (Pp^2)) (λ^4))/(3 (mh^4) (H-Pp) (H+Pp) (ω^2)))

KBtu=((64 (e^2) (mc^2) (mf^2) π ((H^2)+(Pp^2)) (λ^4))/((mh^4) (H-Pp) (H+Pp) (ω^2)))

(((2 (H^2) - (mf^2))/(4(H^2)) (KBtu/. {Pp->(H - (mf^2)/(2H)) })) + O[mf]^4) /. {H-> mc, P-> 0}// Simplify

KFt=-((32 (e^2) (H^2) mc mf π (zlm1-zlm2-zlp1+zlp2) (λ^4))/(Kt (mh^4) Pp))

KFu=-((32 (e^2) π (3 (H^4)+(P^2) (Pp^2)) (zlm1-zlm2-zlp1+zlp2) (λ^4))/(3 Kt (mh^4) Pp))

KFtu=-((32 (e^2) (mc^2) π ((H^2)+(Pp^2)) (zlm1-zlm2-zlp1+zlp2) (λ^4))/(Kt (mh^4) Pp))

KBDirac=((64 (e^2) (mf^2) π (3 (H^4)+(P^2) (Pp^2)) (λ^4))/(3 (mh^4) (H-Pp) (H+Pp) (ω^2)))

KFDirac=-((32 (e^2) (H^2) mc mf π (zlm1-zlm2-zlp1+zlp2) (λ^4))/(Kt (mh^4) Pp))

Self energy part, Collinear log (only from B1, F1)

B1tlog=((32 (e^2) H (3 (H^2)+(P^2)) π (xlm-xlp) (λ^4))/(3 (mh^4) Pp))

B1ulog=((32 (e^2) H (3 (H^2)+(P^2)) π (xlm-xlp) (λ^4))/(3 (mh^4) Pp))

B1tulog=((64 (e^2) H (mc^2) π (xlm-xlp) (λ^4))/((mh^4) Pp))

totBlog=((64 (e^2) H (3 (H^2)-3 (mc^2)+(P^2)) π (xlm-xlp) (λ^4))/(3 (mh^4) Pp))

F1tlog=(1/(3 Kt (mh^4) Pp))16 (e^2) π (3 (H^3) ot (zlm1-zlm2+zlp1-zlp2)+H ot (P^2) (zlm1-zlm2+zlp1-zlp2)+(mf^2) (P^2) (-zlm1+zlm2+zlp1-zlp2)) (λ^4)

F1ulog=(1/(3 Kt (mh^4) Pp))16 (e^2) π (3 (H^3) ot (zlm1-zlm2+zlp1-zlp2)+H ot (P^2) (zlm1-zlm2+zlp1-zlp2)+(mf^2) (P^2) (-zlm1+zlm2+zlp1-zlp2)) (λ^4)

F1tulog=((16 (e^2) (mc^2) π (2 H ot (zlm1-zlm2+zlp1-zlp2)+(mf^2) (-zlm1+zlm2+zlp1-zlp2)) (λ^4))/(Kt (mh^4) Pp))

Simplifying and combining to get  finite contributions (from G + (K_tot_finite) + Bself + Fself).

First, the finite contributions from the G thermal photon and thermal fermion terms:

nB = 1/(Exp[beta ω] -1)

nF = 1/(Exp[beta ot] +1)

G1 = G1t + G1u - G1tu // Simplify

G2 = G2t + G2u - G2tu // Simplify

G3 = G3t + G3u - G3tu // Simplify

G4 = G4t + G4u - G4tu // Simplify

G5 = G5t + G5u - G5tu // Simplify

Ggamma = G1 + G2 + G3 + G4 + G5 /.Pp-> Sqrt[(H^2)-(mf^2)]// Simplify

(1/(15 (mh^8)))128 (e^2) (360 (H^6)-60 (mc^4) (mf^2)+5 (mh^4) (P^2)+8 (mf^2) (P^4)-30 (H^4) (16 (mc^2)+(mh^2)-16 (P^2))-15 (mc^2) ((mh^4)-4 (mf^2) (P^2))+(H^2) ( 120 (mc^4)+15 (mh^4)+30 (mh^2) (P^2)+10 (mc^2) (18 (mf^2)+3 (mh^2)-8 (P^2))+24 (P^2) (-5 (mf^2)+(P^2)))) π (λ^4)

xcv = (360 (H^6)-60 (mc^4) (mf^2)+8 (mf^2) (P^4)-30 (H^4) (16 (mc^2)-16 (P^2))-15 (mc^2) (-4 (mf^2) (P^2))+(H^2) (120 (mc^4)+10 (mc^2) (18 (mf^2)-8 (P^2))+24 (P^2) (-5 (mf^2)+(P^2)))) /. P-> Sqrt[(H^2)- (mc^2)] // Simplify

GgammaDirac = G1t + G2t + G3t + G4t + G5t /.Pp-> Sqrt[(H^2)-(mf^2)]// Simplify (* t channel only *)

G1f = G1ft + G1fu - G1ftu // Simplify

G2f = G2ft + G2fu - G2ftu // Simplify

G2a = G2at + G2au - G2atu // Simplify

(G2f + G2a) + O[ot]^2 // Simplify (* Both G2 from photon and fermion vanish *)

G3a = G3at + G3au - G3atu // Simplify

Gf = (G1f + G2f + G2a + G3a) + O[ot]^2 // Simplify

GfDirac = (G1ft + G2ft + G2at + G3at) + O[ot]^2 // Simplify (* t channel only *)

Next, the finite contributions from the finite part of K:

Kfinite = K4t + K4u - K4tu // Simplify (* All others have pure soft IR or coll div *)

KfiniteDirac = K4t // Simplify (* All others have pure soft IR or coll div *)

Final omega/ot integration to get cross section: first from G and (K_finite):

TotGgamma = Integrate[Ggamma ω nB, {ω, 0, Infinity}] /. {Pp-> Sqrt[(H^2)-(mf^2)]} // Simplify(* Finite G *)

TotGf = Integrate[Normal@Gf ot nF, {ot, 0, Infinity}]/. {Pp-> Sqrt[(H^2)-(mf^2)]} //  Simplify (* Finite G-fermion *)

TotKfinite = Integrate[Kfinite ω nB, {ω, 0, Infinity}]

Now self-energy finite contributions:

TotGself = Normal@totBfinite /. {Pp-> Sqrt[(H^2)-(mf^2)]} // Simplify

TotGfself = Normal@totFfinite  /. {Pp-> Sqrt[(H^2)-(mf^2)]} // Simplify

Total Majorana finite contribution from G, Gf, Kfinite, Bfinite and Ffinite: (all added with + sign)

TotKfinite

Totfinite1 = TotGgamma + TotGf + TotKfinite + 2 TotGself + 2 TotGfself  // Simplify (* for fermion and antifermion self *)

Totfinite = (Coefficient[(mh^8)Normal@Totfinite1, mh, 4] / ((mh^(4 ))) // Simplify)+ (Coefficient[(mh^8)Normal@Totfinite1, mh, 2] / ((mh^(6 ))) // Simplify)+ (Coefficient[(mh^8)Normal@Totfinite1, mh, 0] / ((mh^(8 ))) // Simplify) (* times factorNLOVirtual *)

Totfiniteswave1 = Totfinite /. {H-> mc, P-> 0}// Simplify

Totfiniteswave = (* Now appears to have a 1/mphi^4 contribution which is negative !! *)  (Coefficient[(mh^8)Normal@Totfiniteswave1, mh, 4] / ((mh^(4 ))) // Simplify)+ (Coefficient[(mh^8)Normal@Totfiniteswave1, mh, 2] / ((mh^(6 ))) // Simplify)+ (Coefficient[(mh^8)Normal@Totfiniteswave1, mh, 0] / ((mh^(8 ))) // Simplify)/. {H-> mc, P-> 0}

To match the coefficient to be the same as for sigma_Real so that we can add them together: multiply by P’/(Sqrt[s] = (2H^2-mf^2)/(4H^2)

TotfinitetomatchNLOReal = (Coefficient[(mh^8)Normal@Totfinite1  (((2(H^2)-(mf^2)))/(4(H^2))), mh, 4] / ((mh^(4 ))) // Simplify)+ (Coefficient[(mh^8)Normal@Totfinite1  (((2(H^2)-(mf^2)))/(4(H^2))), mh, 2] / ((mh^(6 ))) // Simplify)+ (Coefficient[(mh^8)Normal@Totfinite1  (((2(H^2)-(mf^2)))/(4(H^2))), mh, 0] / ((mh^(8 ))) // Simplify)

TotfinitetomatchNLORealswave1 = TotfinitetomatchNLOReal /. {H-> mc, P-> 0}// Simplify

TotfinitetomatchNLORealswave =(Coefficient[(mh^8)Normal@TotfinitetomatchNLORealswave1 + O[mf]^4, mh, 4] / ((mh^(4 ))) // Simplify)+ (Coefficient[(mh^8)Normal@TotfinitetomatchNLORealswave1 + O[mf]^4, mh, 2] / ((mh^(6 ))) // Simplify)+ (Coefficient[(mh^8)Normal@TotfinitetomatchNLORealswave1 + O[mf]^4, mh, 0] / ((mh^(8 ))) // Simplify)/. {H-> mc, P-> 0}

Now for the Dirac case:

TotGgammaDirac = Integrate[GgammaDirac ω nB, {ω, 0, Infinity}]

TotGfDirac = Integrate[Normal@GfDirac ot nF, {ot, 0, Infinity}]

TotKfiniteDirac = Integrate[KfiniteDirac ω nB, {ω, 0, Infinity}]

Total finite part:

TotfiniteDirac =(TotGgammaDirac + TotGfDirac + TotKfiniteDirac + 2 totBDirac + 2totFDirac) /. Pp-> Sqrt[(H^2)-(mf^2)]  // Simplify (* for fermion and antifermion self energy *)

TotfiniteDiracswave1 =TotfiniteDirac /. {H-> mc, P-> 0}// Simplify (* Negative for 1/mphi^4 !*)

TotfiniteDiracswave = (* Now appears to have a 1/mphi^4 contribution !! *)  (Coefficient[(mh^8)Normal@TotfiniteDiracswave1, mh, 4] / ((mh^(4 ))) // Simplify)+ (Coefficient[(mh^8)Normal@TotfiniteDiracswave1, mh, 2] / ((mh^(6 ))) // Simplify)+ (Coefficient[(mh^8)Normal@TotfiniteDiracswave1, mh, 0] / ((mh^(8 ))) // Simplify)/. {H-> mc, P-> 0}

TotfiniteDiractomatchNLOReal = (Coefficient[(mh^8)Normal@TotfiniteDirac  (((2(H^2)-(mf^2)) )/(4(H^2))) + O[mf]^4, mh, 4] / ((mh^(4 ))) // Simplify)+ (Coefficient[(mh^8)Normal@TotfiniteDirac (((2(H^2)-(mf^2)))/(4(H^2))) + O[mf]^4, mh, 2] / ((mh^(6 ))) // Simplify)+ (Coefficient[(mh^8)Normal@TotfiniteDirac  (((2(H^2)-(mf^2)))/(4(H^2))) + O[mf]^4, mh, 0] / ((mh^(8 ))) // Simplify)

TotfiniteDiractomatchNLORealswave1 = (TotfiniteDiractomatchNLOReal /.{H-> mc, P-> 0}) + O[mf]^4 // Simplify

TotfiniteDiractomatchNLORealswave = (Coefficient[(mh^8)Normal@TotfiniteDiractomatchNLORealswave1 , mh, 4] / ((mh^(4 ))) // Simplify)+ (Coefficient[(mh^8)Normal@TotfiniteDiractomatchNLORealswave1, mh, 2] / ((mh^(6 ))) // Simplify)+ (Coefficient[(mh^8)Normal@TotfiniteDiractomatchNLORealswave1, mh, 0] / ((mh^(8 ))) // Simplify)

Collinear divergence. Moved it to a separate file.