The finite cross section for dark matter annihilation chi chi -> f fbar (gamma), summed over virtual and real photon (absorption and emission) diagrams . These are sums of non log terms (from virtual and real) and finite remainder from log terms (from virtual and real) . The soft IR divergences cancel between virtual and real non - log K and Ktilde terms while soft plus collinear divergences cancel between virtual and real log terms . There are three contributions : non - log virtual (including self energy diagrams), non - log real, and log (virtual + real) sums that leave a finite remainder. The self-energy contributions are a part of the virtual set; the soft IR divergent pieces are cancelled and the finite remainder of K and Ktilde labelled as the “real” Ktilde_fin terms.

Combined real and virtual contributions  for final cross section  for dark matter annihilation to order T^2 + T^4 at NLO (order alpha).

I. First, the non-log terms:
    . Virtual 
    . Real 

II. Next, the finite part from real and virtual Collinear  log terms 

III. Sums of above

Various factors and LO cross section (just for the record)

factorLO=(Pp/(128 π (H^2) P)) (* 1/(64 (pi^2) s) (Pp/P) (2pi); s = 4(H^2); 2pi from dphi integration *)

IntLO = (( (λ^(4 )))/(((mh))^4)) ( 8 (H^2)((H^2)- (mc^2)) + (mf^2)(5 (mc^2)- 2(H^2)))

sigmaLO = (1/((12 π s))) (Pp/P) IntLO

IntLODirac = (( (λ^(4 )))/(((mh))^4)) (  (H^2)(4 (H^2)- (mc^2)) - (mf^2)( (H^2) - (mc^2)))

factorNLOVirtual=(Pp/((128 (H^2)) P (16 (π^4))))  (* = Pp /Sqrt[s] factorNLOreal *)
(* (1/ (32 s) (Pp/P) (1/(2pi)^4 coming from
1. 1/(64 (pi^2) s) Pp/P (1/(2pi)^4)x [(2pi) from dphi] x[ (2pi) from 2pi delta n_B or n_F] x
[(1/2) from delta(k^2) -> 1/|2 k0|] ; s = 4(H^2) *)

factorNLOreal=Simplify[(1/((32 P) (2H)(((2 π))^4)))]  (* [1/Flux] [1/8(2pi)^4] which comes from 
1. [1/(2(2pi)^6) from two 1/(2pi)^4 times (2pi) from the phase space of the momentum that is integrated out using the en-mom delta fn, and 1/((2pi)^3 2E )],  
2. then a [2pi from dphi integration of nonthermal particle], 
3. [1/2 from 1/(2 |vec{p'}. vec{k}|) of costheta' delta fn,
4. [(2pi) from 2 pi delta(k^2) of thermal particle
5. [(1/2) from 1/|2 k0| of delta(k^2), with Flux = 4 Sqrt[s] Q *)

factorNLOvirtualg=((Pp ω)/(2 H)) (* factor  by which virtual thermal photon has been multiplied *)

factorNLOvirtualf=((Pp Kt)/(2 H)) (* factor  by which virtual thermal fermion has been multiplied *)

Virtual case :gamma and f,fbar: matchreal indicates that the coefficients of the virtual results are matched to the real ones.

intvirtualtotmatchreal=(1/(9 (beta^2) (H^2) (mh^8)))16 (e^2) (-3 (mf^2) (mh^4) (P^2)+8 (H^6) (3 (mf^2)+16 (P^2))-3 (H^2) (mf^2) ((mh^4)+24 (P^4))+(H^4) (-208 (mf^2) (P^2)+64 (P^4))) (π^3) (λ^4)

intvirtualfatotmatchreal=(1/(45 (beta^4) (H^4) (mh^6)))8 (e^2) (5 (beta^2) (H^2) (9 (mf^2) (mh^2) (P^2)+4 (H^4) (3 (mf^2)+8 (P^2))-3 (H^2) (8 (mh^2) (P^2)+(mf^2) ((mh^2)+12 (P^2)))) (π^3)+42 (-2 (H^2)+(mf^2)) (mh^2) (P^2) (π^5)) (λ^4)

intvirtualgammafamatchreal= (* gamma + fa virtual total *) (1/(45 (beta^4) (H^2) (mh^8)))32 (e^2) (P^2) (π^3) (10 (beta^2) (32 (H^6)-3 (H^2) (mh^4)+4 (H^4) ((mh^2)+4 (P^2)))-21 (mh^4) (π^2)) (λ^4)+(1/(45 (beta^4) (H^4) (mh^8)))(8 (e^2) (5 (beta^2) (H^2) (48 (H^6)+3 (mh^4) (P^2)+4 (H^4) (3 (mh^2)-104 (P^2))-9 (H^2) ((mh^4)+4 (mh^2) (P^2)+16 (P^4))) (π^3)+42 (mh^4) (P^2) (π^5)) (λ^4)) (mf^2)+((O[mf])^3)

intvirtualgammafamatchrealswave = intvirtualgammafamatchreal /.{H-> mc, P-> 0} // Simplify

intvirtualtotDiracmatchreal=(1/(9 (beta^2) (mh^8)))16 (e^2) (24 (H^6)+(H^4) (-36 (mf^2)+64 (P^2))-3 (mf^2) ((mh^4)+4 (P^4))+(H^2) (-80 (mf^2) (P^2)+8 (P^4))) (π^3) (λ^4)

intvirtualfatotDiracmatchreal=(1/(45 (beta^4) (H^4) (mh^6)))2 (e^2) (10 (beta^2) (H^2) (24 (H^6)+3 (mf^2) (mh^2) (P^2)-2 (H^4) (6 (mf^2)+9 (mh^2)-4 (P^2))+3 (H^2) (-2 (mh^2) (P^2)+(mf^2) ((mh^2)-4 (P^2)))) (π^3)-21 (2 (H^2)-(mf^2)) (mh^2) (3 (H^2)+(P^2)) (π^5)) (λ^4)

intvirtualgammafaDiracmatchreal=(1/(45 (beta^4) (H^2) (mh^8)))4 (e^2) (10 (beta^2) (H^2) (48 (H^6)-3 (mh^4) (P^2)+4 (H^4) (3 (mh^2)+32 (P^2))+(H^2) (-9 (mh^4)+4 (mh^2) (P^2)+16 (P^4))) (π^3)-21 (mh^4) (3 (H^2)+(P^2)) (π^5)) (λ^4)-(1/(45 (beta^4) (H^4) (mh^8)))(2 (e^2) (10 (beta^2) (H^2) (144 (H^6)-3 (mh^4) (P^2)+4 (H^4) (3 (mh^2)+80 (P^2))+3 (H^2) (3 (mh^4)+4 (mh^2) (P^2)+16 (P^4))) (π^3)-21 (mh^4) (3 (H^2)+(P^2)) (π^5)) (λ^4)) (mf^2)+((O[mf])^3)

intvirtualgammafaDiracmatchrealswave = intvirtualgammafaDiracmatchreal /.{H-> mc, P-> 0} // Simplify

Real case :gamma and f,fbar:

intrealtot=(1/(135 (beta^4) (H^4) (mh^8)))64 (e^2) (π^3) (15 (beta^2) ((H^2) (mf^2) (mh^4) (P^2)+(H^8) (3 (mf^2)+28 (P^2))+(H^6) ((mf^2) (3 (mh^2)-38 (P^2))+4 (P^2) ((mh^2)-(P^2)))+(H^4) (P^2) (-2 (mh^4)+(mf^2) (-7 (mh^2)+3 (P^2))))+2 (3 (H^8)+3 (H^2) (mf^2) (mh^2) (P^2)+3 (mf^2) (mh^4) (P^2)+(H^4) (P^2) ((mf^2)-(P^2))+(H^6) (-3 (mf^2)+38 (P^2))) (π^2)) (λ^4)

intrealtotswave= intrealtot/.{H-> mc, P-> 0} // Simplify

intrealfatot=-(1/(135 (beta^4) (H^4) (mh^8)))2 (e^2) (40 (beta^2) (H^2) (6 (H^6) (mf^2)+18 (mf^2) (mh^4) (P^2)+4 (H^4) (12 (mh^2) (P^2)+(mf^2) (9 (mh^2)-17 (P^2)))-3 (H^2) (16 (mh^4) (P^2)+(mf^2) (3 (mh^4)+20 (mh^2) (P^2)-2 (P^4)))) (π^3)+7 (48 (H^8)-3 (mf^2) (mh^4) (P^2)+16 (H^6) (9 (mf^2)+2 (P^2))-6 (H^4) (24 (mh^2) (P^2)-8 (P^4)+(mf^2) (9 (mh^2)+8 (P^2)))-9 (H^2) (8 (mh^4) (P^2)+(mf^2) (3 (mh^4)-22 (mh^2) (P^2)))) (π^5)) (λ^4)

intrealfatotswave= intrealfatot/.{H-> mc, P-> 0} // Simplify

intrealgammafa=(1/(135 (beta^4) (H^2) (mh^8)))16 (e^2) (π^3) (120 (beta^2) (H^2) (P^2) (14 (H^4)+(mh^4)-2 (H^2) (P^2))+(-18 (H^6)+276 (H^4) (P^2)+63 (mh^4) (P^2)+2 (H^2) (63 (mh^2) (P^2)-25 (P^4))) (π^2)) (λ^4)+(1/(135 (beta^4) (H^4) (mh^8)))(2 (e^2) (π^3) (40 (beta^2) (30 (H^8)-388 (H^6) (P^2)-6 (H^2) (mh^4) (P^2)+3 (H^4) (3 (mh^4)-8 (mh^2) (P^2)+10 (P^4)))+(-1200 (H^6)+213 (mh^4) (P^2)+(H^4) (378 (mh^2)+400 (P^2))+3 (H^2) (63 (mh^4)-398 (mh^2) (P^2))) (π^2)) (λ^4)) (mf^2)+((O[mf])^3)

intrealgammafaswave =  intrealgammafa/.{H-> mc, P-> 0} // Simplify

intrealtotDirac=(1/(135 (beta^4) (H^4) (mh^8)))8 (e^2) (π^3) (15 (beta^2) (H^2) (24 (H^8)+(mf^2) (mh^4) (P^2)+(H^6) (-36 (mf^2)+24 (mh^2)+64 (P^2))-2 (H^4) (3 (mh^4)+4 (mh^2) (P^2)-4 (P^4)+(mf^2) (6 (mh^2)+40 (P^2)))+(H^2) (-2 (mh^4) (P^2)+(mf^2) (3 (mh^4)-4 (mh^2) (P^2)-12 (P^4))))+2 (84 (H^8)+3 (mf^2) (mh^4) (P^2)+4 (H^4) (P^2) ((mf^2)-(P^2))+3 (H^2) (mf^2) (mh^2) (3 (mh^2)+4 (P^2))+(H^6) (-12 (mf^2)+80 (P^2))) (π^2)) (λ^4)

intrealtotDiracswave =  intrealtotDirac/. {H-> mc, P-> 0} // Simplify

intrealfatotDirac=-(1/(135 (beta^4) (H^4) (mh^8)))2 (e^2) (5 (beta^2) (H^2) (144 (H^6) (mh^2)+27 (mf^2) (mh^4) (P^2)-16 (H^4) (9 (mh^4)+11 (mf^2) (P^2)-3 (mh^2) (P^2))+(H^2) (-48 (mh^4) (P^2)+(mf^2) (9 (mh^4)-96 (mh^2) (P^2)-48 (P^4)))) (π^3)+7 (72 (H^8)+3 (mf^2) (mh^4) (P^2)+4 (H^6) (18 (mf^2)-9 (mh^2)-8 (P^2))+3 (H^4) (-9 (mh^4)-12 (mh^2) (P^2)+8 (P^4)+(mf^2) (6 (mh^2)-8 (P^2)))-9 (H^2) ((mh^4) (P^2)+2 (mf^2) ((mh^4)-3 (mh^2) (P^2)))) (π^5)) (λ^4)

intrealfatotDiracswave = intrealfatotDirac/. {H-> mc, P-> 0} // Simplify

intrealgammafaDirac=(1/(135 (beta^4) (H^2) (mh^8)))2 (e^2) (π^3) (120 (beta^2) (12 (H^8)+(H^2) (mh^4) (P^2)+(H^6) (6 (mh^2)+32 (P^2))+(H^4) (3 (mh^4)-6 (mh^2) (P^2)+4 (P^4)))+(168 (H^6)+63 (mh^4) (P^2)+36 (H^4) (7 (mh^2)+24 (P^2))+(H^2) (189 (mh^4)+252 (mh^2) (P^2)-200 (P^4))) (π^2)) (λ^4)-(1/(135 (beta^4) (H^4) (mh^8)))(2 (e^2) (π^3) (5 (beta^2) (432 (H^8)+15 (H^2) (mh^4) (P^2)+16 (H^6) (9 (mh^2)+49 (P^2))-3 (H^4) (9 (mh^4)+16 (mh^2) (P^2)-32 (P^4)))+(600 (H^6)-3 (mh^4) (P^2)+2 (H^4) (63 (mh^2)-100 (P^2))-6 (H^2) (33 (mh^4)-47 (mh^2) (P^2))) (π^2)) (λ^4)) (mf^2)+((O[mf])^3)

intrealgammafaDiracswave= intrealgammafaDirac/.{H-> mc, P-> 0} // Simplify

Collinear contributions, sum over real and virtual: gamma and f,fbar:

intCollgammaRV=(1/(135 (beta^4) (H^4) (mh^6)))16 (e^2) (15 (beta^2) (H^2) (-((mf^2) (mh^2) (P^2))+12 (H^4) ((mf^2)+4 (P^2))+(H^2) (8 (mh^2) (P^2)+(mf^2) (3 (mh^2)-56 (P^2)))) (π^3)+2 (-((mf^2) (mh^2) (P^2))+3 (H^4) ((mf^2)-8 (P^2))-(H^2) (8 (mh^2) (P^2)+(mf^2) (3 (mh^2)+23 (P^2)))) (π^5)) (λ^4)

intCollgammaRVswave = intCollgammaRV/.{H-> mc, P-> 0} // Simplify

intCollfaRV=-(1/(135 (beta^4) (H^4) (mh^6)))2 (e^2) (60 (beta^2) (H^2) (6 (H^4) (mf^2)+5 (mf^2) (mh^2) (P^2)-(H^2) (24 (mh^2) (P^2)+(mf^2) (9 (mh^2)+38 (P^2)))) (π^3)+7 (-357 (mf^2) (mh^2) (P^2)+6 (H^4) (5 (mf^2)+24 (P^2))-(H^2) (168 (mh^2) (P^2)+(mf^2) (63 (mh^2)+58 (P^2)))) (π^5)) (λ^4)

intCollfaRVswave= intCollfaRV/.{H-> mc, P-> 0} // Simplify

intCollgammafaRV=(1/(135 (beta^4) (H^4) (mh^6)))2 (e^2) (π^3) (60 (beta^2) (-7 (H^2) (mf^2) (mh^2) (P^2)+6 (H^6) (3 (mf^2)+16 (P^2))+(H^4) (40 (mh^2) (P^2)+(mf^2) (15 (mh^2)-74 (P^2))))+(2483 (mf^2) (mh^2) (P^2)-6 (H^4) (27 (mf^2)+232 (P^2))+(H^2) (1048 (mh^2) (P^2)+(mf^2) (393 (mh^2)+38 (P^2)))) (π^2)) (λ^4)

intCollgammafaRVswave =  intCollgammafaRV/. {H-> mc, P-> 0} // Simplify

intCollgammaRVDirac=-(1/(135 (beta^4) (H^4) (mh^6)))8 (e^2) (π^3) (-15 (beta^2) (24 (H^8)-(H^2) (mf^2) (mh^2) (P^2)+6 (H^6) (-3 (mf^2)+(mh^2)+4 (P^2))+(H^4) (2 (mh^2) (P^2)+(mf^2) (3 (mh^2)-26 (P^2))))+(24 (H^6)-(mf^2) (mh^2) (P^2)+12 (H^4) ((mf^2)+(mh^2)+2 (P^2))+(H^2) (4 (mh^2) (P^2)+(mf^2) (9 (mh^2)+28 (P^2)))) (π^2)) (λ^4)

intCollgammaRVDiracswave= intCollgammaRVDirac/.{H-> mc, P-> 0} // Simplify

intCollfaRVDirac=(1/(270 (beta^4) (H^4) (mh^6)))(e^2) (60 (beta^2) (H^2) (12 (H^4) (2 (mf^2)+3 (mh^2))-7 (mf^2) (mh^2) (P^2)+(H^2) (12 (mh^2) (P^2)+5 (mf^2) (3 (mh^2)+8 (P^2)))) (π^3)+7 (-144 (H^6)+147 (mf^2) (mh^2) (P^2)+12 (H^4) ((mf^2)+21 (mh^2)-12 (P^2))+(H^2) (84 (mh^2) (P^2)+11 (mf^2) (63 (mh^2)+4 (P^2)))) (π^5)) (λ^4)

intCollfaRVDiracswave= intCollfaRVDirac/. {H-> mc, P-> 0} // Simplify

intCollgammafaRVDirac=(1/(270 (beta^4) (H^4) (mh^6)))(e^2) (60 (beta^2) (H^2) (96 (H^6)-11 (mf^2) (mh^2) (P^2)+(H^4) (-48 (mf^2)+60 (mh^2)+96 (P^2))+(H^2) (20 (mh^2) (P^2)+(mf^2) (27 (mh^2)-64 (P^2)))) (π^3)+(-1392 (H^6)+1045 (mf^2) (mh^2) (P^2)-12 (H^4) (9 (mf^2)-131 (mh^2)+116 (P^2))+(H^2) (524 (mh^2) (P^2)+(mf^2) (4707 (mh^2)-140 (P^2)))) (π^5)) (λ^4)

intCollgammafaRVDiracswave = intCollgammafaRVDirac/.{H-> mc, P->0} // Simplify

Total = "virtual" + "real" + "Collinear real+virtual" :

inttotal=(1/(135 (beta^4) (H^2) (mh^8)))32 (e^2) (π^3) (120 (beta^2) (H^2) (P^2) (15 (H^4)+(mh^4)+(H^2) (4 (mh^2)+3 (P^2)))+(-9 (H^6)+138 (H^4) (P^2)+34 (mh^4) (P^2)-(H^2) (24 (mh^2) (P^2)+25 (P^4))) (π^2)) (λ^4)+(1/(135 (beta^4) (H^4) (mh^8)))(4 (e^2) (20 (beta^2) (H^2) (102 (H^6)-12 (mh^4) (P^2)+(H^4) (45 (mh^2)-1012 (P^2))+3 (H^2) (6 (mh^4)-63 (mh^2) (P^2)-62 (P^4))) (π^3)+(-600 (H^6)+1600 (mh^4) (P^2)+4 (H^4) (27 (mh^2)+50 (P^2))+(H^2) (291 (mh^4)-578 (mh^2) (P^2))) (π^5)) (λ^4)) (mf^2)+((O[mf])^3)

inttotalswave= inttotal/. {H-> mc, P-> 0} // Simplify

inttotalDirac=(1/(135 (beta^4) (H^2) (mh^8)))8 (e^2) (60 (beta^2) (H^2) (18 (H^6)+(mh^4) (P^2)+12 (H^4) ((mh^2)+4 (P^2))+(H^2) (3 (mh^4)+4 (mh^2) (P^2)+6 (P^4))) (π^3)+(42 (H^6)+17 (mh^4) (P^2)-24 (H^4) ((mh^2)-9 (P^2))+(H^2) (51 (mh^4)-24 (mh^2) (P^2)-50 (P^4))) (π^5)) (λ^4)-(((e^2) (40 (beta^2) (H^2) (648 (H^6)+15 (mh^4) (P^2)+4 (H^4) (45 (mh^2)+338 (P^2))-3 (H^2) (9 (mh^4)-36 (mh^2) (P^2)-64 (P^4))) (π^3)+(2400 (H^6)-1309 (mh^4) (P^2)+(H^4) (612 (mh^2)-800 (P^2))+(H^2) (-6255 (mh^4)+1268 (mh^2) (P^2))) (π^5)) (λ^4)) (mf^2))/(270 ((beta^4) (H^4) (mh^8)))+((O[mf])^3)

inttotalDiracswave= inttotalDirac/.{H-> mc, P-> 0} // Simplify