Version X,10.10.24, with corrected l for u-channel diagrams

<<FeynCalc`

01ttFS

Theory and Necessary Information

(1) Model for the process in hand is as follows

ℒ =  - (1/4) (F^(μν ))(Subscript[F, μν])+ (OverBar[f])(((iγ)^μ)(Subscript[D, μ]) - (Subscript[m, f]))f  + (1/2) (OverBar[χ]) (((iγ)^μ)(∂_(μ ))- (Subscript[m, χ]))χ + ((Subscript[D, μ])ϕ )^†((D^μ)ϕ) - ((Subscript[m, ϕ])^2)(SuperDagger[ϕ])ϕ + ( λ (OverBar[χ])(Subscript[P, L])(SuperMinus[f])(SuperPlus[ϕ])  +  (SuperDagger[λ])(SuperPlus[(OverBar[f])])(Subscript[P, R])χ (SuperMinus[ϕ]))

Definitions and Approximations

Primary definitions

definition of b((Subscript[p, i]),(Subscript[p, f]),k) = b(x,y,z)

b[x_,y_,z_] :=(1/2) (((SP[2x-z,2y-z])/((
SP[z]-2*SP[x,z])*((SP[z]-2*SP[y,z]))))+((SP[2x+z,2y+z])/((SP[z]+2*SP[x,z])*((SP[z]+2*SP[y,z])))))

(Subscript[Q, f])= -e;(Subscript[Q, f']) = -e;(Subscript[Q, ϕ])=-e; (* If all charges are taken to be the same, then f+ and phi+ imply that for the t (u)-channel, scalar momentum flows upwards (downwards), so that its t-channel momeuntum is l = q-p and u-channel is l' = q'-p; I had taken q-p' earlier, which is wrong. Then for antiparticle with momentum p, the Feynman rule is taken as momentum going backwards, i.e.,as -p. So the overall signs will be wrong in Diagrams 1_u, 3_u and possibly 4_u *)

Rotated from CM frame in which q’ lies along z to one where p’ along z by rotating by -theta
Component form of  4-vectors q’, q, p’, p and k

q’	=   	(	H, 	-P Sin θ , 0,	 P Cos θ)
q	=   	(	H, 	P Sin θ, 0, 	-P Cos θ )
p’	=   	(	H’,  0, 0, P’ )
p	=   	(	H’, 0, 0, -P’)
k	=   	(	ω, K Sin θ’ Cos ϕ’, K Sin θ’ Sin ϕ’, K Cos θ’ )

Expansion of various scalar products (Notations : Ct = Cos θ, Ctp = Cos θ’

Scalar products

SP[k]=0;
SP[p] = (mf^2);
SP[pp] = (mf^2);
SP[q] = (mc^2);
SP[qp] = (mc^2);
SP[l] = SP[q-p];(* dot products of these not yet defined *)
SP[lp] = SP[qp-p]; (* dot products of these not yet defined *)
SP[t] = (mt^2); (* for thermal fermion; otherwise thermal fermion divergences; corrected this compared to the earlier version *)

(3.3)Matrix elements

(3.3.1)t- channel LO Matrix element

* SP[q-p]-((Subscript[m, ϕ])^2)   -> lphi ; l = q - p

Mt = (I*λ)(I*ComplexConjugate[λ])(SpinorUBar[pp,mf].GA[6].SpinorU[qp,mc])(SpinorVBar[q,mc].GA[7].SpinorV[p,mf])((I/lphi)) // ExpandScalarProduct

(3.3.2) t- channel NLO Matrix element

Note : only 11 part of the whole propagator for all of the relevant particles are taken. 12,21 and 22 parts for the same are not included.
* Following matrix element is written without δ((k^2)) and action of δ((k^2)) is taken by taking (k^2)->0. 
* It contains only G-polarization part of Integrand inside t- channel NLO matrix element.
* SP[q-p]-((Subscript[m, ϕ])^2)   -> -((Subscript[m, ϕ])^2)	and	SP[q-p+k]-((Subscript[m, ϕ])^2)   -> -((Subscript[m, ϕ])^2); for fermion thermal, replaced dummy t by dummy k for convenience

MGtot8Ft =(I*ComplexConjugate[λ])(I*λ)(-I*(Subscript[Q, f]))(-I*(Subscript[Q, f]))(SpinorUBar[pp,mf].GA[μ].SpinorV[p,mf])(SpinorVBar[q,mc].GA[7].(-DiracSlash[qp+q-t]+mf).GA[ν].(DiracSlash[t]+mf).GA[6].SpinorU[qp,mc])((I/lkphi))(((-I)/((4(H^2)- 2 SP[q+qp,t]))))((I/(4(H^2))))(-MT[μ,ν])  (* Priyotosh's diagram, with thermal fermion; actually k = t - qp *)

MGtot8At = (I*ComplexConjugate[λ])(I*λ)(-I*(Subscript[Q, f]))(-I*(Subscript[Q, f]))(SpinorUBar[pp,mf].GA[μ].SpinorV[p,mf])(SpinorVBar[q,mc].GA[7].(-DiracSlash[t]+mf).GA[ν].(DiracSlash[qp+q-t]+mf).GA[6].SpinorU[qp,mc])((I/lkphi))(* actually k = q-t *) (((-I)/((4(H^2)- 2 SP[q+qp,t]))))((I/((4(H^2)))))(-MT[μ,ν]) (* Priyotosh's diagram, with thermal anti-fermion *)

MGtot9Ft =(I*ComplexConjugate[λ])(I*λ)(-I*(Subscript[Q, f]))(-I*(Subscript[Q, f]))(SpinorUBar[pp,mf].GA[μ].SpinorV[p,mf])(SpinorVBar[q,mc].GA[7].(-DiracSlash[t]+mf).GA[6]. SpinorU[qp,mc])(-FV[qp-q+2t,ν])((I/lkphi))((I/lkphi))(((-I)/(4(H^2))))(-MT[μ,ν])  (* Priyotosh's diagram modified with photon emitted from scalar instead of fermion, with thermal fermion; lkphi may not be correct, but using heavy scalar approx, so doesn't matter *)

(3.3.1)u- channel LO Matrix element

* SP[q'-p]-((Subscript[m, ϕ])^2)   -> lphip ; l’ = q’ - p

Mu =(I*ComplexConjugate[λ])(I*λ)(SpinorUBar[pp,mf].GA[6].SpinorU[q,mc]) (SpinorVBar[qp,mc].GA[7].SpinorV[p,mf])((I/lphip)) // ExpandScalarProduct

u- channel NLO Matrix element

Note : only 11 part of the whole propagator for all of the relevant particles are taken. 12,21 and 22 parts for the same are not included.
* Following matrix element is written without δ((k^2)) and action of δ((k^2)) is taken by taking (k^2)->0. 
* It contains only G-polarization part of Integrand inside t- channel NLO matrix element.
* SP[q'-p]-((Subscript[m, ϕ])^2)   -> lphip -> -((Subscript[m, ϕ])^2) ; Note q’-p p’-q for virtual only, not real case.

MGtot8Fu =(I*ComplexConjugate[λ])(I*λ)(-I*(Subscript[Q, f]))(-I*(Subscript[Q, f]))(SpinorUBar[pp,mf].GA[μ].SpinorV[p,mf])(SpinorVBar[qp,mc].GA[7].(-DiracSlash[qp+q-t]+mf).GA[ν].(DiracSlash[t]+mf).GA[6].SpinorU[q,mc])((I/lkphi))(((-I)/((4(H^2)- 2 SP[q+qp,t]))))((I/(4(H^2))))(-MT[μ,ν])  (* Priyotosh's diagram, with thermal fermion; actually k = t - qp *)

MGtot8Au = (I*ComplexConjugate[λ])(I*λ)(-I*(Subscript[Q, f]))(-I*(Subscript[Q, f]))(SpinorUBar[pp,mf].GA[μ].SpinorV[p,mf])(SpinorVBar[qp,mc].GA[7].(-DiracSlash[t]+mf).GA[ν].(DiracSlash[qp+q-t]+mf).GA[6].SpinorU[q,mc])((I/lkphi))(* actually k = q-t *) (((-I)/((4(H^2)- 2 SP[q+qp,t]))))((I/((4(H^2)))))(-MT[μ,ν]) (* Priyotosh's diagram, with thermal anti-fermion *)

MGtot9Fu =(I*ComplexConjugate[λ])(I*λ)(-I*(Subscript[Q, f]))(-I*(Subscript[Q, f]))(SpinorUBar[pp,mf].GA[μ].SpinorV[p,mf])(SpinorVBar[qp,mc].GA[7].(DiracSlash[t]+mf).GA[6]. SpinorU[q,mc])(FV[qp-q+2t,ν])((I/lkphi))((I/lkphi))(((-I)/(4(H^2))))(-MT[μ,ν])  (* Priyotosh's diagram modified with photon emitted from scalar instead of fermion, with thermal fermion; lkphi may not be correct, but using heavy scalar approx, so doesn't matter *)

(3.4) Component form of  4-vectors q’, q, p’, p and k

q’	=   	(	H, 	-P.Sin(θ) ,0,	 P. Cos(θ)	)

q	=   	(	H, 	P.Sin(θ), 	0, 	-P Cos(θ)	)

p’	=   	(	H’,  	0, 		0, 	P’		)

p	=   	(	H’,  	0, 		0, 	-P’		)

k	=   	(ω,  	K.Cos(ϕ’).Sin(θ’),   K.Sin(θ’).Sin(ϕ’),   K.Cos(θ’))

(3.5) expansion of various scalar products (Notations : Ct = Cos(θ), Ctp = Cos(θ’) FRAME ROTATED; run this AFTER traces

SP[p,q]=H*H-P*Pp*Ct;
SP[p,qp]=H*H+Pp*P*Ct;
SP[p,pp]=H*H+Pp*Pp;

SP[q,qp]=H*H+P*P;
SP[q,pp]=H*H+P*Pp*Ct;

SP[qp,pp]=H*H-P*Pp*Ct;

SP[k,p]=ω*H+Pp*K*Ctp;
SP[k,q] = ω*H-K*P*St*Stp*Cos[phip]+K*P*Ct*Ctp;
SP[k,qp] = ω*H+K*P*St *Stp *Cos[phip]-K*P*Ct*Ctp;
SP[k,pp]= ω*H- K*Pp*Ctp;

SP[t,p]=ot*H+Pp*Kt*Ctp; (* ot^2 = Kt^2 + mf^2 *)
SP[t,q] = ot*H-Kt*P*St*Stp*Cos[phip]+Kt*P*Ct*Ctp;
SP[t,qp] = ot*H+Kt*P*St *Stp *Cos[phip]-Kt*P*Ct*Ctp;
SP[t,pp]= ot*H- Kt*Pp*Ctp;

(* Non-dynamical scalar*)
(* For the virtual case *)
lphi = -(mh^2) (* lphi = l^2 - mphi^2 = (q-p)^2 - mphi^2 *);
lphip = -(mh^2) (* lphip = lp^2 - mphi^2 = (qp-p)^2 - mphi^2 *);
lkphi = -(mh^2) (* lkphi = (l+k)^2 - mphi^2 *);
lkphip = -(mh^2) (* lkphip = (lp+k)^2 - mphi^2 *);
(* For the real case *)
(* lphi = -(mphi^2); lphi = l^2 - mphi^2 = (q-p)^2 - mphi^2 *)
(* lphip = -(mphi^2); lphip = lp^2 - mphi^2 = (qp-p)^2 - mphi^2 *)
lmkphi = -(mh^2) (* lmkphi = (l-k)^2 - mphi^2 *);
lmkphip = -(mh^2) (* lmkphip = (lp-k)^2 - mphi^2 *);

Trace Calculations

(4.1) t channel contributions

Mo00P = (Mt) *ComplexConjugate[Mt]+(Mu)*ComplexConjugate[Mu] - (Mt) *ComplexConjugate[Mu]-(Mu)*ComplexConjugate[Mt]//Contract//FermionSpinSum//DiracSimplify//Simplify

Mt2 =(Mt) *ComplexConjugate[Mt]//Contract//FermionSpinSum//DiracSimplify//Simplify

MGtoto8Ft1 = (MGtot8Ft) *ComplexConjugate[Mt]+(Mt)*ComplexConjugate[MGtot8Ft]//Contract//FermionSpinSum//DiracSimplify//Simplify

MGtoto8Ft = MGtoto8Ft1 /. mt-> mf  // Simplify

MGtoto8At1 = (MGtot8At) *ComplexConjugate[Mt]+(Mt)*ComplexConjugate[MGtot8At]//Contract//FermionSpinSum//DiracSimplify//Simplify

MGtoto8At = MGtoto8At1 /. mt-> mf  // Simplify

MGtoto9Ft1 = (MGtot9Ft) *ComplexConjugate[Mt]+(Mt)*ComplexConjugate[MGtot9Ft]//Contract//FermionSpinSum//DiracSimplify//Simplify

ComplexExpand[Re[MGtoto9Ft1]]

MGtoto9Ft = ComplexExpand[Re[MGtoto9Ft1]]/.  mt-> mf  // Simplify (* There are not enough directions in 2-> 2 process for the Levi Civita tensors to survive *)

u-channel contributions:

MGtoto8Fu1 = (MGtot8Fu) *ComplexConjugate[Mu]+(Mu)*ComplexConjugate[MGtot8Fu]//Contract//FermionSpinSum//DiracSimplify//Simplify

MGtoto8Fu = MGtoto8Fu1 /.   mt-> mf  // Simplify

MGtoto8Au1 = (MGtot8Au) *ComplexConjugate[Mu]+(Mu)*ComplexConjugate[MGtot8Au]//Contract//FermionSpinSum//DiracSimplify//Simplify

MGtoto8Au = MGtoto8Au1 /.   mt-> mf  // Simplify

MGtoto9Fu1 = (MGtot9Fu) *ComplexConjugate[Mu]+(Mu)*ComplexConjugate[MGtot9Fu]//Contract//FermionSpinSum//DiracSimplify//Simplify

MGtoto9Fu =  ComplexExpand[Re[MGtoto9Fu1]]/.  mt-> mf  // Simplify

Cross terms in tu

MGtoto8Ftu1 = (MGtot8Ft) *ComplexConjugate[Mu]+(Mu)*ComplexConjugate[MGtot8Ft]+(MGtot8Fu) *ComplexConjugate[Mt]+(Mt)*ComplexConjugate[MGtot8Fu]//Contract//FermionSpinSum//DiracSimplify//Simplify (* tu and ut both *)

MGtoto8Ftu = MGtoto8Ftu1 /.  mt-> mf  // Simplify

MGtoto8Atu1 = (MGtot8At) *ComplexConjugate[Mu]+(Mu)*ComplexConjugate[MGtot8At]+(MGtot8Au) *ComplexConjugate[Mt]+(Mt)*ComplexConjugate[MGtot8Au]//Contract//FermionSpinSum//DiracSimplify//Simplify

MGtoto8Atu = MGtoto8Atu1 /.  mt-> mf  // Simplify

MGtoto9Ftu1 = (MGtot9Ft) *ComplexConjugate[Mu]+(Mu)*ComplexConjugate[MGtot9Ft]+(MGtot9Fu) *ComplexConjugate[Mt]+(Mt)*ComplexConjugate[MGtot9Fu]//Contract//FermionSpinSum//DiracSimplify//Simplify

MGtoto9Ftu = ComplexExpand[Re[MGtoto9Ftu1]] /. mt-> mf  // Simplify

Now for the t channel F- terms:

MGtoto8Gt = MGtoto8Ft /.  ot -> -ot// Simplify

MGtoto8gt = MGtoto8At /.  ot -> -ot// Simplify

MGtoto9Gt = MGtoto9Ft /.  ot -> -ot// Simplify

Now the u-channel F- terms:

MGtoto8Gu = MGtoto8Fu /.  ot -> -ot// Simplify

MGtoto8gu = MGtoto8Au /.  ot -> -ot// Simplify

MGtoto9Gu = MGtoto9Fu /.  ot -> -ot// Simplify

The tu cross terms for F-

MGtoto8Gtu = MGtoto8Ftu /.  ot -> -ot // Simplify

MGtoto8gtu = MGtoto8Atu /.  ot -> -ot// Simplify

MGtoto9Gtu = MGtoto9Ftu /.  ot -> -ot// Simplify

(4.4) Replacing K by (k^0)  ((k^(0 )) = ω)  for thermal photons only

t-channel terms:

FGtoto81t = MGtoto8Ft //Simplify

FGtoto82t = MGtoto8Gt //Simplify

AGtoto81t = MGtoto8At //Simplify

AGtoto82t = MGtoto8gt //Simplify

FGtoto91t = MGtoto9Ft //Simplify

FGtoto92t = MGtoto9Gt //Simplify

The u-channel terms:

FGtoto81u = MGtoto8Fu //Simplify

FGtoto82u = MGtoto8Gu //Simplify

AGtoto81u = MGtoto8Au //Simplify

AGtoto82u = MGtoto8gu //Simplify

FGtoto91u = MGtoto9Fu //Simplify

FGtoto92u = MGtoto9Gu //Simplify

The tu cross terms:

FGtoto81tu = MGtoto8Ftu //Simplify

FGtoto82tu = MGtoto8Gtu //Simplify

AGtoto81tu = MGtoto8Atu //Simplify

AGtoto82tu = MGtoto8gtu //Simplify

FGtoto91tu = MGtoto9Ftu //Simplify

FGtoto92tu = MGtoto9Gtu //Simplify

(5) Phase Space Integral Calculations with special choice of frame of reference for minimizing calculations

(5.1) ϕ’ integration

IntG8tf1ttot = Integrate[Integrate[(FGtoto81t/.Stp-> Sqrt[1-Ctp Ctp]),{phip, 0,2*π}],{Ctp, -1, 1}] // Simplify (* FGtoto81t = MGtoto8Ft //Simplify   *)

IntG8tf2ttot = Integrate[Integrate[(FGtoto82t/.Stp-> Sqrt[1-Ctp Ctp]),{phip, 0,2*π}],{Ctp, -1, 1}] // Simplify (* FGtoto81t = MGtoto8Ft //Simplify   *)

IntG8ta1ttot = Integrate[Integrate[(AGtoto81t/.Stp-> Sqrt[1-Ctp Ctp]),{phip, 0,2*π}],{Ctp, -1, 1}] // Simplify

IntG8ta2ttot = Integrate[Integrate[(AGtoto82t/.Stp-> Sqrt[1-Ctp Ctp]),{phip, 0,2*π}],{Ctp, -1, 1}] // Simplify

IntG9tf1ttot = Integrate[Integrate[(FGtoto91t/.Stp-> Sqrt[1-Ctp Ctp]),{phip, 0,2*π}],{Ctp, -1, 1}] // Simplify (* FGtoto81t = MGtoto8Ft //Simplify   *)

IntG9tf2ttot = Integrate[Integrate[(FGtoto92t/.Stp-> Sqrt[1-Ctp Ctp]),{phip, 0,2*π}],{Ctp, -1, 1}] // Simplify (* FGtoto81t = MGtoto8Ft //Simplify   *)

IntG8tf1utot = Integrate[Integrate[(FGtoto81u/.Stp-> Sqrt[1-Ctp Ctp]),{phip, 0,2*π}],{Ctp, -1, 1}] // Simplify (* FGtoto81t = MGtoto8Ft //Simplify   *)

IntG8tf2utot = Integrate[Integrate[(FGtoto82u/.Stp-> Sqrt[1-Ctp Ctp]),{phip, 0,2*π}],{Ctp, -1, 1}] // Simplify (* FGtoto81t = MGtoto8Ft //Simplify   *)

IntG8ta1utot = Integrate[Integrate[(AGtoto81u/.Stp-> Sqrt[1-Ctp Ctp]),{phip, 0,2*π}],{Ctp, -1, 1}] // Simplify

IntG8ta2utot = Integrate[Integrate[(AGtoto82u/.Stp-> Sqrt[1-Ctp Ctp]),{phip, 0,2*π}],{Ctp, -1, 1}] // Simplify

IntG9tf1utot = Integrate[Integrate[(FGtoto91u/.Stp-> Sqrt[1-Ctp Ctp]),{phip, 0,2*π}],{Ctp, -1, 1}] // Simplify (* FGtoto81t = MGtoto8Ft //Simplify   *)

IntG9tf2utot = Integrate[Integrate[(FGtoto92u/.Stp-> Sqrt[1-Ctp Ctp]),{phip, 0,2*π}],{Ctp, -1, 1}] // Simplify (* FGtoto81t = MGtoto8Ft //Simplify   *)

IntG8tf1tutot = Integrate[Integrate[(FGtoto81tu/.Stp-> Sqrt[1-Ctp Ctp]),{phip, 0,2*π}],{Ctp, -1, 1}] // Simplify (* FGtoto81t = MGtoto8Ft //Simplify   *)

IntG8tf2tutot = Integrate[Integrate[(FGtoto82tu/.Stp-> Sqrt[1-Ctp Ctp]),{phip, 0,2*π}],{Ctp, -1, 1}] // Simplify (* FGtoto81t = MGtoto8Ft //Simplify   *)

IntG8ta1tutot = Integrate[Integrate[(AGtoto81tu/.Stp-> Sqrt[1-Ctp Ctp]),{phip, 0,2*π}],{Ctp, -1, 1}] // Simplify

IntG8ta2tutot = Integrate[Integrate[(AGtoto82tu/.Stp-> Sqrt[1-Ctp Ctp]),{phip, 0,2*π}],{Ctp, -1, 1}] // Simplify

IntG9tf1tutot = Integrate[Integrate[(FGtoto91tu/.Stp-> Sqrt[1-Ctp Ctp]),{phip, 0,2*π}],{Ctp, -1, 1}] // Simplify (* FGtoto81t = MGtoto8Ft //Simplify   *)

IntG9tf2tutot = Integrate[Integrate[(FGtoto92tu/.Stp-> Sqrt[1-Ctp Ctp]),{phip, 0,2*π}],{Ctp, -1, 1}] // Simplify (* FGtoto81t = MGtoto8Ft //Simplify   *)

(5.2) θ Integration

Int00P = Integrate[(Mo00P/. St->Sqrt[1-Ct Ct]),{Ct, -1, 1}] // Simplify

IntG8f1ttot =  Integrate[(IntG8tf1ttot/. St->Sqrt[1-Ct Ct]),{Ct, -1, 1}] // Simplify

IntG8f2ttot =  Integrate[(IntG8tf2ttot/. St->Sqrt[1-Ct Ct]),{Ct, -1, 1}] // Simplify

IntG8tf = (IntG8f1ttot + IntG8f2ttot)/. {Pp-> Sqrt[(H^2)-(mf^2)],Kt-> Sqrt[(ot^2)-(mf^2)]} // Simplify

IntG8a1ttot =  Integrate[(IntG8ta1ttot/. St->Sqrt[1-Ct Ct]),{Ct, -1, 1}] // Simplify

IntG8a2ttot =  Integrate[(IntG8ta2ttot/. St->Sqrt[1-Ct Ct]),{Ct, -1, 1}] // Simplify

IntG8ta = (IntG8a1ttot + IntG8a2ttot)/. {Pp-> Sqrt[(H^2)-(mf^2)],Kt-> Sqrt[(ot^2)-(mf^2)]} // Simplify

IntG8tswave = ((IntG8tf + IntG8ta)/. {H-> mc, P-> 0})  // Simplify

IntG9f1ttot =  Integrate[(IntG9tf1ttot/. St->Sqrt[1-Ct Ct]),{Ct, -1, 1}] // Simplify

IntG9f2ttot =  Integrate[(IntG9tf2ttot/. St->Sqrt[1-Ct Ct]),{Ct, -1, 1}] // Simplify

IntG9tf = (IntG9f1ttot + IntG9f2ttot)/. {Pp-> Sqrt[(H^2)-(mf^2)],Kt-> Sqrt[(ot^2)-(mf^2)],(LC[p, pp, qp,t]-> 0)} // Simplify

IntG9tfswave = IntG9tf/. {H-> mc, P-> 0}  // Simplify

IntG8f1utot =  Integrate[(IntG8tf1utot/. St->Sqrt[1-Ct Ct]),{Ct, -1, 1}] // Simplify

IntG8f2utot =  Integrate[(IntG8tf2utot/. St->Sqrt[1-Ct Ct]),{Ct, -1, 1}] // Simplify

IntG8uf = (IntG8f1utot + IntG8f2utot)/. {Pp-> Sqrt[(H^2)-(mf^2)],Kt-> Sqrt[(ot^2)-(mf^2)]} // Simplify

IntG8a1utot =  Integrate[(IntG8ta1utot/. St->Sqrt[1-Ct Ct]),{Ct, -1, 1}] // Simplify

IntG8a2utot =  Integrate[(IntG8ta2utot/. St->Sqrt[1-Ct Ct]),{Ct, -1, 1}] // Simplify

IntG8ua = (IntG8a1utot + IntG8a2utot)/. {Pp-> Sqrt[(H^2)-(mf^2)],Kt-> Sqrt[(ot^2)-(mf^2)]} // Simplify

IntG8uswave = ((IntG8uf + IntG8ua)/. {H-> mc, P-> 0})  // Simplify

IntG9f1utot =  Integrate[(IntG9tf1utot/. St->Sqrt[1-Ct Ct]),{Ct, -1, 1}] // Simplify

IntG9f2utot =  Integrate[(IntG9tf2utot/. St->Sqrt[1-Ct Ct]),{Ct, -1, 1}] // Simplify

IntG9uf = (IntG9f1utot + IntG9f2utot)/. {Pp-> Sqrt[(H^2)-(mf^2)],Kt-> Sqrt[(ot^2)-(mf^2)],((OverBar[ϵ])^((OverBar[p]) (OverBar[pp]) (OverBar[qp]) (OverBar[t]))) -> Eps} // Simplify

IntG9ufswave = ((IntG9uf )/. {H-> mc, P-> 0})  // Simplify

IntG8f1tutot =  Integrate[(IntG8tf1tutot/. St->Sqrt[1-Ct Ct]),{Ct, -1, 1}] // Simplify

IntG8f2tutot =  Integrate[(IntG8tf2tutot/. St->Sqrt[1-Ct Ct]),{Ct, -1, 1}] // Simplify

IntG8tuf = (IntG8f1tutot + IntG8f2tutot)/. {Pp-> Sqrt[(H^2)-(mf^2)],Kt-> Sqrt[(ot^2)-(mf^2)]} // Simplify

IntG8a1tutot =  Integrate[(IntG8ta1tutot/. St->Sqrt[1-Ct Ct]),{Ct, -1, 1}] // Simplify

IntG8a2tutot =  Integrate[(IntG8ta2tutot/. St->Sqrt[1-Ct Ct]),{Ct, -1, 1}] // Simplify

IntG8tua = (IntG8a1tutot + IntG8a2tutot)/. {Pp-> Sqrt[(H^2)-(mf^2)],Kt-> Sqrt[(ot^2)-(mf^2)]} // Simplify

IntG8tuswave = ((IntG8tuf + IntG8tua)/. {H-> mc, P-> 0}) // Simplify

IntG8tswave - IntG8tuswave // Simplify

IntG9f1tutot =  Integrate[(IntG9tf1tutot/. St->Sqrt[1-Ct Ct]),{Ct, -1, 1}] // Simplify

IntG9f2tutot =  Integrate[(IntG9tf2tutot/. St->Sqrt[1-Ct Ct]),{Ct, -1, 1}] // Simplify

IntG9tuf = (IntG9f1tutot + IntG9f2tutot)/. {Pp-> Sqrt[(H^2)-(mf^2)],Kt-> Sqrt[(ot^2)-(mf^2)]} // Simplify

IntG9tufswave = ((IntG9tuf )/. {H-> mc, P-> 0})  // Simplify

IntG8swave = IntG8tswave + IntG8uswave - IntG8tuswave // Simplify

IntG8 = ((IntG8tf + IntG8ta) +(IntG8uf + IntG8ua)-(IntG8tuf + IntG8tua) )/. mc-> Sqrt[(H^2)-(P^2)]// Simplify

IntG8swave = IntG8 /.  {H-> mc, P-> 0}  // Simplify

IntG8matchreal = Normal@(Normal@((IntG8 )(2(H^2) - (mf^2))/(4(H^2))) + O[mf]^3) + O[ot]^4 // Simplify (* x Pp / root(s); T2 term is helicity suppressed; all terms vanish in s-wave limit  *)

IntG8swavematchreal = (IntG8matchreal /.{H-> mc, P-> 0})  // Simplify

IntG8swaveT2matchreal  = Normal@(IntG8swavematchreal)// Simplify

IntG8tmatchreal = Normal@(Normal@((IntG8tf + IntG8ta )(2(H^2) - (mf^2))/(4(H^2))) + O[mf]^3) + O[ot]^4 // Simplify (* x Pp / root(s) *)

IntG8tmatchrealswave =  IntG8tmatchreal /. {H-> mc, P-> 0} // Simplify

(IntG9tf)/. mc-> Sqrt[(H^2)-(P^2)] // Simplify

IntG9  = (IntG9tf + IntG9uf - IntG9tuf)/. mc-> Sqrt[(H^2)-(P^2)] // Simplify

IntG9matchreal = Normal@(Normal@((IntG9 )(2(H^2) - (mf^2))/(4(H^2))) + O[mf]^3) + O[ot]^4 // Simplify (* x Pp / root(s); T2 term is helicity suppressed *)

IntG9tmatchreal = Normal@(Normal@((IntG9tf )(2(H^2) - (mf^2))/(4(H^2))) + O[mf]^3) + O[ot]^4 // Simplify (* x Pp / root(s) *)

IntG9tmatchrealswave = IntG9tmatchreal /. {H-> mc, P-> 0} // Simplify

IntG9swavematchreal = IntG9matchreal /. {H-> mc, P-> 0} // Simplify

Integrating over ot : Int...matchreal means that the virtual has the same overall coeff as real terms (called Int’ in paper)

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

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

I2 = Integrate [(ω nB), {ω, 0, Infinity}]

I4 = Integrate [((ω^3) nB), {ω, 0, Infinity}]

If2 = Integrate [(Kt nF/. Kt-> ot), {ot, 0, Infinity}]

If4 = Integrate [(Kt (ot^2)nF/. Kt-> ot), {ot, 0, Infinity}]

IntG8Intmatchreal = Coefficient[IntG8matchreal, ot, 0] If2 + Coefficient[IntG8matchreal, ot, 2] If4 // Simplify

IntG8Intmatchrealswave = IntG8Intmatchreal/. {H-> mc, P-> 0} // Simplify

IntG8DiracIntmatchreal = Coefficient[IntG8tmatchreal, ot, 0] If2 + Coefficient[IntG8tmatchreal, ot, 2] If4 // Simplify

IntG8DiracIntmatchrealswave =  IntG8DiracIntmatchreal/. {H-> mc, P-> 0} // Simplify

IntG9Intmatchreal = Coefficient[IntG9matchreal, ot, 0] If2 + Coefficient[IntG9matchreal, ot, 2] If4 // Simplify

IntG9Intmatchrealswave = IntG9Intmatchreal/. {H-> mc, P-> 0} // Simplify

IntG9DiracIntmatchreal = Coefficient[IntG9tmatchreal, ot, 0] If2 + Coefficient[IntG9tmatchreal, ot, 2] If4 // Simplify

IntG9DiracIntmatchrealswave =  IntG9DiracIntmatchreal/. {H-> mc, P-> 0} // Simplify