NLO thermal corrections to dark matter annihilation cross sections: a novel approach
Prabhat Butola D. Indumathi Pritam Sen




1,2]Prabhat Butola[prabhatb@imsc.res.in,
0009-0003-2824-455X]

1,2]D. Indumathi[indu@imsc.res.in,
0000-0001-6685-4760]

3]Pritam Sen[pritam.sen@tifr.res.in,
0000-0003-0751-5560]



[1]Homi Bhabha National Institute, Mumbai, India

[2]Institute of Mathematical Sciences, Chennai, India

[3]Tata Institute of Fundamental Research, Mumbai, India

NLO thermal corrections to dark matter
annihilation cross sections: a novel approach
    [
    April 29, 2024
====================================================================================




The dark matter relic density has been increasingly accurately measured by
successive generations of experiments. The Boltzmann equation determines
the yields using the dark matter annihilation cross section as one of
the inputs; the accurate computation of the latter including thermal
contributions thus assumes importance. We report here the next-to-leading
order (NLO) thermal corrections to the cross sections for (Majorana)
dark matter annihilation to standard model fermions: χχ→ f f, via charged scalars. We use a novel approach,
utilising the technique of Grammer and Yennie, extended to thermal field
theories, where the cancellation of soft infra-red divergences occurs
naturally. We present the NLO thermal cross sections in full detail for
both the relativistic case as well as in the non-relativistic limit. Our
independent calculation verifies earlier results where the leading
contribution at order O(T^2) was shown to be proportional to
the square of the fermion mass in the non-relativistic limit, just as
at leading order. We find that the O(T^4) contributions have
the same dependence on the fermion mass as well.




Keywords: Thermal field theory, Dark Matter annihilation, NLO cross
section



§ INTRODUCTION



Evidence for the presence of dark matter (DM) in our Universe arises
from several observations, primarily due to its
interaction with standard model particles via the gravitational channel
<cit.>. Evidence for the existence of dark
matter includes its effect in gravitational

lensing, structure formation and its presence is a necessity for
stability of rotating galaxies; various experiments have set limits
or constrained various properties of dark matter <cit.> via whole-sky or galactic observations.

Cosmic microwave background radiation is sensitive to the baryonic,
dark and the total matter densities respectively, through the height of
the peaks in its power spectrum. The relic abundance of dark matter (DM)
is encoded in its relic density, currently measured to be Ω_c h^2
= 0.1200 ± 0.0012 by the PLANCK collaboration <cit.>,
where h is the reduced Hubble constant, h = H_0/100. This is higher
than the abundance of baryonic matter and lower than the share of dark
energy in the energy density of our Universe.

Unfortunately, there is as yet no direct evidence for dark matter via,
for example, whole-sky searches of self-annihilation of dark matter
to produce gamma rays <cit.>, galactic surveys
<cit.>, etc. Various non-baryonic and baryonic
DM candidates have been postulated <cit.> such as WIMPs,
axion-like particles (ALP), neutralino dark matter, sterile neutrinos,
particles in extra dimensional theories, primordial black holes, etc.;
see Refs. <cit.> for reviews. A
large sub-set of theoretical proposals and experiments to detect
them from cosmology, colliders, etc., exist; a few recent ones can
be found in Refs. <cit.>.

One of the popular theories of dark matter consists of a SUSY-inspired model of bino-like cold dark matter
produced thermally in the Universe post-inflation. The number density
of the dark matter is then determined by Boltzmann equations
dependent on the interaction rate (annihilation/production of the dark
matter) as well as the expansion rate of the Universe. The dark matter
is in thermal, chemical and kinetic equilibrium with standard model (SM)
particles until its interaction rate falls below the Hubble expansion
rate <cit.>. In such a freeze-out mechanism for thermal
production of dark matter <cit.>, the DM starts decoupling from the background plasma,
when the Hubble rate of expansion of our Universe became comparable to
the annihilation rate of DM species via, for example, χχ→ f
f. The temperature where DM decouples from SM particles
(f), depends on the mass of DM as well as the thermally averaged
annihilation cross section ⟨σ v_Mø l⟩, and we
have a range of temperature where the mass of DM can be comparable
to decoupling temperature, i.e., x ≡ m_χ/T ∼O(1). In this range of temperature, thermal corrections can
play an important role in determining the precise relic density of dark
matter. In particular, the larger the dark matter annihilation cross
section, the later the DM goes out of equilibrium, leading to a lower
relic density, and vice versa.

A contrasting model is when the dark matter is never in equilibrium with
the SM particles; instead, the coupling of DM to SM particles is so
small that the annihilation of χχ→ f f
can be neglected. The amount of dark matter in the Universe then keeps
on increasing due to the reverse annihilation process f f→χχ until the dark matter freezes-in to the present relic
density <cit.>. 

In any case, present-day determination of dark matter relic densities
are becoming so precise that these cross sections must necessarily
be calculated to next-to-leading order (NLO). While NLO calculations
of various relevant cross sections exist <cit.>, these mostly
focussed on higher order quantum corrections to the leading cross
sections in various models. It was only recently that the higher order
thermal corrections to these cross sections have been calculated
<cit.>. In particular,
Beneke et al. <cit.> computed the thermal NLO corrections
to the dark matter annihilation cross section in a bino-like model of
dark matter using thermal field theory that we describe briefly below. In
this model, the dark matter candidate is an SU(2) × U(1) singlet
Majorana fermion χ which interacts with SM doublet fermions, f =
(f^0,f^-)^T, via scalar partners ϕ = (ϕ^+ , ϕ^0)^T through a
Yukawa interaction, 
    L    = -1/4 F_μν
    F^μν + f( i D - m_f ) f +
    1/2χ( i ∂ - m_χ)
    χ
            + (D^μϕ)^†(D_μϕ) -
      m_ϕ^2 ϕ^†ϕ + ( λχ P_L f^-
      ϕ^+ +  h.c.) .

The thermal corrections to the annihilation of dark matter via
χχ→ f f at NLO in this model were studied in
Ref. <cit.> where the cancellation of the infrared
divergences was explicitly demonstrated. Subsequently, an all-order
proof of the cancellation of infra-red divergences in such thermal field
theories containing scalars and fermions interacting with photons in a
heat bath at finite temperature was given in Refs. <cit.>, using a generalised Grammer and Yennie technique
<cit.> although the finite
remainder was not calculated.  The advantage of this approach is the
automatic cancellation of the soft infra-red (IR) divergences so that we
can directly compute the finite remainder. The calculations are simpler
because they are explicitly IR finite.  Here we compute the IR-finite
contributions to the annihilation cross section at NLO, using the same
technique.  We calculate both the O(T^2) and O(T^4)
corrections both in the relativistic case and in the non-relativistic
limit using this novel approach to obtain an independent calculation of
the NLO thermal contributions to the annihilation cross section, and
determine the contribution in the non-relativistic limit. The latter
are important at freeze-out when m_χ/T ∼ 20. We show that the
O(T^2) terms are proportional to the square of the fermion
masses, just as for the LO cross section, a result which arises due
to helicity suppression of the process due to the presence of Majorana
fermions. This result was first obtained in Ref. <cit.>
who subsequently calculated <cit.> the cross section
for annihilation of Dirac-type dark matter particles using the Operator
Product Expansion approach.

We make a careful study of the dependence of both the O(T^2)
and O(T^4) corrections; the O(T^2) terms are indeed
proportional to m_f^2 in the non-relativistic limit, while there
are terms proportional to the square of the dark matter momenta p^2
in the general case, which of course are small in the non-relativistic
limit. We find in addition that the leading 
O(T^4) terms are also proportional to the square of the
fermion masses in the non-relativistic limit; hence it appears that
helicity conservation determines this behaviour beyond the leading order
as well. We discuss the general behaviour of the cross section as well
as the behaviour in the non-relativistic limit; detailed expressions
for the cross section are available online <cit.> as sets of
Mathematica Notebooks.

In the next section, we highlight some key points of thermal field
theory in the real time formalism. In Section 3, we present the
details of the Grammer and Yennie approach which allows us to easily
separate the IR-finite part of the cross section. We also discuss some
technical details about the various contributions which allow us to
simplify the calculation at NLO. In Section 4 we calculate the leading
order cross section, which is presented for completeness. In Sections 5
and 6, we calculate the NLO cross section in two different scenarios.
For both the LO and NLO calculations, we use the FeynCalc
<cit.> software
with Mathematica 13 <cit.>. We end with Discussions in
Section 7.



§ THERMAL FIELD THEORY



We briefly review the real-time formulation of thermal field theories
<cit.>.
The ensemble average of an operator can be written <cit.> as
the expectation value of time-ordered products in the thermal vacuum. The
unique feature is that the integration in the complex time plane is
defined over a special path shown in Fig. <ref>, from an
initial time, t_i to a final time, t_i - i β, where β is
the inverse temperature of the heat bath, β = 1/T;. As a result,
the thermal fields satisfy the periodic boundary conditions,

    φ(t_0) = ±φ(t_0 - i β) ,

where the sign +1 (-1) corresponds to boson (fermion) fields.
The real-time axis C_1 allows for scattering to be defined while the
remaining portions of the contour lead to a doubling of the fields,
which can be of type-1 (the physical fields), or type-2 (the “ghost”
fields), that live on C_1 and C_2 respectively.



Only type-1 fields can occur on external legs <cit.> while
fields of both types can occur on internal legs, with the off-diagonal
elements of the propagator allowing for conversion of one type into
another. The Feynman rules (for the propagators and vertices) for the
Lagrangian density given in Eq.<ref> in a heat bath at finite
temperature T are listed in Appendix <ref>. We merely note
here that all propagators (which are of 2 × 2 matrix form), can
be written as the sum of two terms, one which is temperature-independent and the other which contains the explicitly thermal
dependence which we call the thermal part.



§ THE GRAMMER AND YENNIE APPROACH



The application of the generalised Grammer and Yennie approach for the
case of the bino-like Lagrangian described in Eq. <ref> has been
described in detail in Refs. <cit.>. Here we
summarise the main results for the sake of completeness, clarity and
use of notation.

We use a generalised version of the Grammer and Yennie technique
<cit.> which was used for the
zero temperature quantum field theory to separate the infrared (IR)
divergences in the theory. The technique was extended to the case of
thermal field theory, first for the case of the interaction of charged
fermions with photons in thermal QED <cit.>, and later
for both charged thermal scalars and a general theory of charged scalars
and fermions interacting with dark matter particles with thermal QED
corrections <cit.>.  In all these instances,
the all-order cancellation of infra-red (IR) divergences was shown by
separating the virtual photon propagator into two parts, the so-called K
and G photons, which we define below, such that the IR divergence
was completely contained in the K photon contribution and the G
photon contribution was IR finite. However, as mentioned earlier, the
IR-finite part was not explicitly calculated in Ref. <cit.>
and is the goal of this work.

The Grammer and Yennie procedure is to insert a virtual photon with
momentum k into a lower order graph, for (say) the process A(p) B(q)
→ C(p'), such that the momentum q, flowing into a given vertex V,
is not soft; this special vertex V is labelled as a hard vertex with
respect to which the graph can be separated into a notional (initial)
p_i-leg and a (final) p_f-leg. We will identify V, p_i and
p_f for the dark matter annihilation diagrams of interest in the next
section. The K and G parts are subsequently identified with respect
to these definitions.  The additional inserted photon can be soft, in
which case the points of insertion will be labelled as soft vertices,
and will contribute to the IR divergence of the corresponding cross
section. The procedure starts with re-arranging, in its propagator term,

    -i
    g^μν   →-i{[ g^μν - b_k (p_f, p_i) k^μ
            k^ν] + [ b_k (p_f, p_i) k^μ k^ν]
            } , 
       ≡-i{[ G_k^μν] +
     [K_k^μν] } ,

where

    b_k(p_f,p_i) = 1/2[ (2p_f-k) ·
    (2p_i-k)/((p_f-k)^2-m^2)((p_i-k)^2-m^2) + (k ↔ -k)
    ] ,

is defined symmetrically in k→ -k for the thermal case (in contrast
to the original definition in Ref. <cit.>), and is a
function of k as well as the momenta, p_f, p_i.
It can then be shown that the K photon insertions contain the IR
divergent pieces while the G photon contribution is IR finite.



 §.§ The K photon insertion and the IR divergent piece


We briefly describe, for the sake of completeness, the separation and
cancellation of IR divergences in the theory. Details are available in
Refs. <cit.>. Let us start with the
example of insertion of an additional virtual photon in the lowest-order
diagrams shown in Fig. <ref>. The key step is to specifically
define a “hard vertex” V with respect to which we can identify the
momenta p_f and p_i. We (arbitrarily) define this special vertex V
to be the one where the momentum q' (of χ) enters. Then, a higher
order diagram can be obtained by inserting a virtual photon anywhere in
either of the LO diagrams shown in Fig. <ref>. For every such
insertion, p_f (p_i) will correspond to the momentum p' or p
depending on whether the final (initial) vertex of the virtual photon
was inserted on the p' or p leg.  Then for example, a virtual photon
inserted with one vertex on the final fermion leg and one on the scalar
leg (as in Diagram (1) of Fig. <ref>), will have a virtual photon
propagator involving the term b_k(p',p) since the p_f line corresponds
to the p' leg, while the p_i line corresponds to the scalar leg
from V to X together with the anti-fermion leg with momentum p.
Thus the corresponding K and G photon insertions can be defined
using Eq. <ref>. For the u-channel diagram, such an insertion
will involve b_k(p',p'); see Ref. <cit.> for more details.



We now demonstrate the factorisation of the IR divergences in the
K-photon insertions. We take the example of the case when the inserted
photon has both its vertices, μ, ν, on the p' fermion line. The
relevant part of the matrix element can be expressed as (note that there
is an integration measure over the photon momentum, k)

    M_NLO^p'p', K photon   ∝
    b_k(p',p') k^μ  k^ν [u(p',m_f) γ_μ S_ fermion^t_μ,t_ν(p'+k, m_f) γ_ν S_ fermion^t_ν,t_V (p', m_f)
    P_R   u(q',m_χ) ] ⋯ , 
    
        = b_k(p',p')  [u_p' k S^μ,ν_p'+k k S^ν,V_p'
    P_R   u_q'] ⋯ ,

where the ellipsis refer to terms independent of k and we have
lightened the notation for convenience in the second line, with the
superscripts on S referring to the thermal type.  The term
in the square brackets can be simplified, using the identities given in
Appendix <ref> as

    [ ]     = [u_p' k {S^μ,ν_p'+k kS^ν,V_p'} 
    P_R   u_q'] , 
    
        = (-1)^t_ν+1 [u_p' k {S_p'^μ,V δ_t_ν,t_V -
    S_p'+k^μ,V δ_t_ν,t_μ}
    P_R   u_q'] , 
    
        = (-1)^t_ν+1 [u_p' { 0 -[(-1)^t_μ+1 δ_t_ν,t_μ] }
    P_R   u_q'] ,

where the first term vanishes since the integrand is odd in k → -k
while both b_k and the measure are even under this exchange. The crucial
step is seen in the second line of Eq. <ref>, where the Feynman
identities (see Eq. <ref>) are used to reduce the portions with
k insertions to a difference of two terms. Hence the
matrix element factorises into two terms, one that is proportional to
the lower order matrix element, viz.,

    M_NLO^p'p', K photon   ∝ b_k(p',p')   M_LO ,

and the other containing the IR divergence.



  §.§.§ Insertion into a general n^ th order graph


Consider an n^ th order graph with n virtual/real photon
vertices; for specificity, we consider higher order corrections to
the t-channel diagram. When an additional virtual K-photon is
inserted into this graph, there are several possible locations where the
additional photon vertices can be inserted. Adding all the contributions
gives differences of two terms as seen above, with, more importantly,
sets of terms cancelling against each other until only one term, that
is proportional to the lower order matrix element, is left.

A similar result is obtained when the photon vertices are inserted so
that one vertex μ is on the p' leg and the other ν vertex is
on the p leg; that is, either on the scalar line or the anti-fermion
line. Here again, it turns out that the K photon contribution is
proportional[Technically, a virtual photon insertion leads to
the number of vertices increasing by 2; we have used the index (n+1)
to indicate that it is an n^ th order graph with an additional
virtual photon.] to the lower order matrix element:

    M_n+1^p'p, K photon   ∝ b_k(p',p)   M_n ,

In fact, such an insertion is the sum of contributions when the
second ν vertex of the virtual photon is inserted on the scalar,
and when it is inserted on the anti-fermion line. These contributions
separately cancel among themselves, leaving one term in the former and
two in the latter. These cancel across the vertex V, leaving behind
a contribution that is again proportional to the lower order graph, as
seen in Eq. <ref> above. In short, it is found <cit.>
that the total matrix element due to insertion of the virtual K photon
into an n^ th order diagram is given by

    M_n+1^K photon, tot =     [ ie^2/2∫ d^4 k/(2π)^4 {0pt12ptδ_t_μ,t_1 δ_t_ν,t_1  D^t_μ,t_ν (k)  [ b_k(p',p') + b_k(p,p) ] . . 
           +. . δ_t_μ,t_V δ_t_ν,t_V 
    	 D^t_μ,t_ν (k)  [ -2b_k(p',p) ] 0pt12pt}0pt16pt]
    	 M_n , 
    ≡    [B ] M_n ,

where the prefactor containing the IR divergence can be expressed as,

    B     = ie^2/2∫ d^4 k/(2π)^4  D^11 (k)  [ b_k(p',p') - 2 b_k(p',p) + b_k(p,p) ] , 
    	 
       ≡ie^2/2∫ d^4 k/(2π)^4 
    	D^11 (k)  [ J^2(k) ] ,

since the thermal types of the hard/external vertices must be type-1. We
see that each term is proportional to the (11) component of the photon
contribution and this is crucial for achieving the cancellation between
virtual and real photon insertions. When the contributions are summed
over all orders in perturbation theory, the IR divergent factor B
exponentiates and cancels a similar IR divergent factor arising from
real soft photon insertions; for details please refer to
Ref. <cit.>.



  §.§.§ Insertion of the photon into a thermal fermion line



We already know that the presence of the bosonic number operator,
n_B(| k^0 |), makes the IR contribution to the cross section
potentially IR divergent while the nature of the fermionic number
operator, n_F(| p^0 |), results in a finite IR contribution to
the cross section. In particular, we shall see below that the K-photon
contribution when the additional photon is inserted into a thermal
fermion (that is, when the second, temperature-dependent term in the
fermionic propagator contributes; see Eq. <ref>), is zero. We
will use this result later in the computation of the NLO cross section
of interest. In order to show this, we begin with the insertion of one
of the vertices of the additional K photon on any thermal fermion line
at the vertex μ, lying between vertices μ_q+1 and μ_q;
see Fig. <ref>.



There are now two fermion propagators, one between vertices μ_q and μ
and the other between vertices μ and μ_q+1. The relevant part of
the matrix element is given by (where we have chosen to insert into the
p' fermion line for specificity),

    M_NLO^ thermal f, Kγ   ∝ b(p_f, p_i) k^μ k^ν[ u_p'γ_μ_1⋯γ_μ_q S^q,μ_p'+∑_qγ_μ S^μ,q+1_p'+∑_q+k γ_μ_q+1] ⋯ ,
    
        = b(p_f, p_i) k^ν[ u_p'γ_μ_1⋯γ_μ_q {S^q,μ_p'+∑_qk S^μ,q+1_p'+∑_q+k}γ_μ_q+1]
    ⋯ ,

where the momentum flowing in the fermion line between vertices q
and (q+1) into which the additional photon was inserted is p' +
l_1 + ⋯ + l_q ≡ p' + Σ_q. The term in curly braces
can be simplified using the property of the fermion propagator shown
in Eq. <ref> where S is also defined:

    {⋯}    = [(p' + Σ_q + m_f) k
    (p' + Σ_q + k + m_f) ]
    S^q,μ_p'+∑_qS^μ,q+1_p'+∑_q+k ,
    
        = [ (2(p' + Σ_q) · k + k^2)  (p'
    + Σ_q + m_f) - ((p'+Σ_q)^2- m_f^2) k] , 
        = 0  for thermal fermions.

In particular, if both the fermion propagators on either side of the
K-photon insertion at the vertex μ are thermal, then both
terms in Eq. <ref> vanish due to the delta function
terms δ ((p'+Σ_q)^2- m_f^2) and δ ((p'+Σ_q+k)^2-
m_f^2) in the thermal part of the respective propagators. We will use
this result in the next section when we compute the NLO cross section.

In order to show the IR divergence cancellation between virtual and
real photon contributions, a similar treatment of the real photon
corrections to the process is done, by separating the polarisation
sums for the insertion of real photons into so-called K
and G contributions:

    [                                           ∑_ polϵ^*,μ (k) ϵ^ν (k)                                                                 =                                                         -
    	g^μν ,;                                                                                                                                   = -{[ g^μν -
    	b̃_k(p_f,p_i) k^μ k^ν] +
    	[b̃_k(p_f,p_i) k^μ k^ν] } ,;                                                                                                                                   ≡                                       -{[ G_k^μν] + [
     K_k^μν]} . ]

Since k^2=0 for real photons, we define,

    b̃_k(p_f, p_i) = b_k(p_f, p_i)|_k^2 = 0 =
    	p_f ·p_i/ k·p_f  k·p_i ,

where p_i (p_f) is the momentum p' or p depending on whether the
real photon insertion was on the p' or p leg in the n^ th
order matrix element
M_n+1 (or its conjugate M_n+1^†). Again,
the insertion of a K real photon into an
n^ th order graph leads to a cross section that is proportional
to the lower order one, with the cancellation again occurring pair-wise,
to give

    |M_n+1^Kγ, tot|^2
    	   ∝ -e^2 [b̃_k(p,p) -2
    	b̃_k(p',p) + b̃_k(p',p') ] ,
    		
       ≡ -e^2 J^2(k) .

The contributions from the virtual K (Eq. <ref>) and real
K (Eq. <ref>) photon corrections are IR
divergent and cancel order by order in the theory. We are therefore left
with the IR-finite G photon contribution, that we will now evaluate
in Sections <ref> and <ref>. However,
first, for completeness, we present the LO results.



§ THE DARK MATTER ANNIHILATION CROSS SECTION AT LO



The leading order (LO) contribution to the annihilation process χχ→ f f arises from the t- and u-channel
processes shown in Fig. <ref>. The momenta in the CM frame are
q', q → p', p, with the choices

    q'^μ    = (H, 0, 0, P) ,        
    q^μ    = (H, 0, 0, -P) , 
    
    p'^ μ    = (H, P'sinθ, 0, P'cosθ) ,        
    p^μ    = (H, -P'sinθ, 0, -P'cosθ) ,

where the centre of momentum energy, s = 4H^2, and θ is the angle
between the initial and final momenta (q⃗', p⃗'). The cross
section at leading order for this 2→ 2 process is given as usual by,

    σ_LO    =    1/64π^2 s |p⃗'⃗|/|q⃗'⃗|∫dΩ[∑_spins |M_LO^t - M_LO^u |^2 ] , 
        =    1/32π s P'/P∫dcosθ[∑_spins |M_LO^t - M_LO^u |^2 ] ,

where the integration over
the azimuth ϕ is trivial, and a summation over both
final state and initial state spins is to be performed since all
helicity configurations contribute to the total cross section.

The contribution from the matrix element can be expressed as

    ∫dcosθ∑_spins |M_LO^t -
    M_LO^u |^2 
       ≡    Int^t_LO + Int^u_LO - Int^tu_LO , 
       ≡    Int_LO ,

where the subscript LO denotes the LO contribution and the three terms
on the right hand side refer to the square of the t-channel matrix
element, the square of the u-channel matrix element, and the tu cross
terms respectively. Note that the cross term vanishes for Dirac-type
dark matter particles. The t-, u-channel matrix elements at LO
are given by,

    M_LO^t     = 
    i λ ^2(v(q,m_χ) P_L   v(p,m_f))  Δ(l)
    
    (u(p',m_f) P_R   u(q',m_χ) ) , 
    M_LO^u     =
    i λ ^2(v(q',m_χ) P_L   v(p,m_f))  Δ(l')
    
    (u(p',m_f) P_R   u(q,m_χ) ) ,
 
where we have lightened the notation for the scalar propagator so that
Δ(l) = i/(l^2-m_ϕ^2), and similarly for l', with
l=q-p ≡ p' - q' and l'=q' -p ≡ p'-q referring to
the momentum of the intermediate scalar for the t and u channels
respectively, and P_R,L≡ (1 ±γ_5)/2. Adding all the
contributions, we find,

    Int_LO    = 2 λ ^4/P P'(2 H^2-m_Φ^2){-[H^2 (m_χ^2-2
    m_Φ^2)+m_χ^2 P'^2+m_Φ^4]
    log(2 H^2-m_Φ^2+2 P P')/(2 H^2-m_Φ^2-2 P P'). 
    
       .     
    +4 P P' (2 H^2-m_Φ^2)
    (2 H^4-2 H^2 m_Φ^2+m_Φ^4-2 P^2 P'^2)/4
    H^4-4 H^2 m_Φ^2+m_Φ^4-4 P^2 P'^2} ,

where m_Φ^2 ≡ m_χ^2 + m_f^2 - m_ϕ^2. The logarithmic
terms arise when either of p, p' are
collinear with one of q, q'; however, these are not divergent;
they contribute at a single phase space point and are well-behaved. Then
the invariant cross section at LO can be written as

    σ_LO(s)     =   1/32 π sP'/P Int_LO ,

where P'^2 = H^2 - m_f^2, P^2 = H^2 - m_χ^2, and 4 H^2 =
s, the usual Mandelstam variable.



  
LO cross section in the non-relativistic limit:
In the non-relativistic limit, when the velocity of the dark matter
particles is small, we can write P = m_χ v and H^2 ≈
m_χ^2(1+v^2). Since P ≪ H, m_Φ, we can expand the log term as

    log(2 H^2-m_Φ^2+2 P P')/2 H^2-m_Φ^2-2 P P'P small⟶ 
    2 [ 2 P P'/2 H^2-m_Φ^2 + 2/3(2 P P'/2 H^2-m_Φ^2)^3 + ⋯] ,

to get

    Int_LO   v small⟶  8 λ^4 m_χ^2 m_f^2/(m_χ^2+ m_ϕ^2 - m_f^2)^2 + O(v^2) ,

Notice that the cross section is proportional to the square of the fermion
mass, which is a well-known result <cit.>.



Hence, the LO cross section can be written as

    σ_LO(s)     =   1/32 π sP'/P Int_LO , 
       v small⟶λ^4/4 π sP'/P[m_χ^2 m_f^2/(m_χ^2+ m_ϕ^2 - m_f^2)^2 + O(v^2)
    ] .

This term is the usual velocity-independent “a” term in s  σ(s)
v_rel = a + b v^2 for the annihilation process <cit.>;
with the relative velocity between the two dark matter particles given
by v_rel = 2 v, where v is the CM velocity of either of the
particles. Terms of order O(v^2) can be calculated by
retaining higher orders in the expansion of Eq. <ref>.
 
We can repeat the calculation in the limit when the scalar is
much heavier than the other particles, viz., m_ϕ^2 ≫ m_χ^2
≳ m_f^2. Then l^2 ≡ (q-p)^2 ≪ m_ϕ^2 (where we have
implicitly assumed that √(s)≪ m_ϕ) and the scalar
propagator can be approximated by i D_ϕ = i/(l^2-m_ϕ^2)
→ -i/m_ϕ^2, so that we get

    Int_LO^ heavy scalar    =
    8 λ^4/3 m_ϕ^4[6 H^4-3 H^2 m_χ^2+P'^2 (2 P^2-3 m_χ^2)] , 
        = 8 λ^4/3 m_ϕ^4[8 H^2( H^2 - m_χ^2) +m_f^2 (5 m_χ^2 - 2 H^2)] ,

where we have substituted for P, P' in the last line. This gives us a
cross section,

    σ_LO^ heavy scalar    =
    1/12π sP'/Pλ^4/m_ϕ^4[8 H^2( H^2 - m_χ^2) +m_f^2 (5 m_χ^2 - 2 H^2)] .

The first term in the square brackets is proportional to H^2 = m_χ^2
≡ P^2 and the
second is proportional to m_f^2. In the non-relativistic limit, with
H^2 = m_χ^2(1+ v^2), P= m_χ v, this matches the expression
given in Eq. <ref>, with the further approximation,
(m_χ^2+m_ϕ^2-m_f^2) → m_ϕ^2, which is valid in the heavy
scalar limit. The important features are instantly visible within this
approximation: the 1/m_ϕ^4 dependence arising from the scalar
propagator, and the m_f^2 dependence of the cross section in the
non-relativistic limit, H→ m_χ. Therefore, we present the NLO
calculation of the cross section in the next sections in two parts: in
the next section, we present the results in the heavy-scalar
limit, where the expressions are shorter and the features can be easily
understood. In the subsequent section, we then present the results with
the fully dynamical scalar propagator. While the detailed results are
cumbersome and can be found on-line <cit.>, the results in the
non-relativistic limit have been presented for comparison and discussion
here.



§ THE DARK MATTER ANNIHILATION CROSS SECTION AT
NLO—“HEAVY-SCALAR” APPROXIMATION



We apply
the Grammer and Yennie technique <cit.>, generalised for
the case of thermal field theories as explained in
Section <ref>. There are two sets of contributions at higher
order: insertions of a virtual photon in the LO diagrams shown in
Fig. <ref>, as well as insertions of real photons which can
be both emitted into or absorbed from the heat bath at temperature
T.  Each of these contributions can be factorised into K and G
(or K and G) parts for virtual (or real)
photon insertions respectively, by applying the Grammer and Yennie
technique. As a consequence, the IR divergences are contained in the K
(K) contributions respectively and furthermore, were shown
earlier to cancel not only at NLO <cit.>, but at all
orders <cit.> in perturbation theory. Note that only soft
photon emissions/absorptions of the process χχ→ f f
(γ) are contained in the K contribution. The hard
real-photon emissions/absorptions that contribute to G
photon insertions do not contribute to the process of interest
here. Hence, by applying the Grammer and Yennie technique, the NLO
corrections to χχ→ f f can simply be calculated
by computing the G terms of the virtual photon contributions alone.

The various contributions[Note that fermion self-energy
corrections do not contribute to the G photon terms.] at NLO are
shown in Fig. <ref>.  In order to compute the G photon
contribution, we will replce the photon propagator by its G photon
part; see Eq. <ref>.  Furthermore, throughout this section, we
will approximate the scalar propagator by i D_ϕ = i/(l^2-m_ϕ^2)
→ -i/m_ϕ^2, and lift this approximation in the next section.
Analogous to the LO cross section, we can write the next-to-leading order
(NLO) virtual contribution as

    σ_NLO   ∝   [∑_spins(M_LO^t - M_LO^u
    )^†(M_NLO^t - M_NLO^u ) + h.c.
    ] ,

where each of the higher order t-channel terms in Eq. <ref>
gets contributions from the five t-channel diagrams shown in
Fig. <ref> (and similarly, their u-channel counterparts).

We
will discuss the details of the calculation using Diagram 1 as an
example, before we list the contributions from all the remaining
diagrams.





 §.§ The NLO contribution from Diagram 1



We start by discussing the pure t-channel NLO contribution from the
first diagram in Fig. <ref>. Since the scalar is very heavy, the
thermal contribution from the scalar propagators, which contain
a factor n_ϕ = 1/(exp[β E_ϕ]-1), can be dropped since
E_ϕ≳ m_ϕ≫ m_χ and β m_χ≡ x ∼ 20
near freeze-out. Hence the scalar propagators can only be of (11) or
(22) type; see Appendix <ref> for details on the thermal
and temperature independent parts of various propagators. Since all
external lines can only be of type-1 <cit.>, and scalar
propagators of (12) or (21) type can be neglected, it can be seen that
all propagators are of type (11) and all vertices are of type-1
alone. In fact, we will find that this holds for all the contributions
at NLO, for all the diagrams.

Therefore, there are three possibilities for Diagram 1, viz.,




  * the photon propagator (with momentum k) is thermal, 


  * the fermion propagator (with momentum p' +k) is thermal,


  * or both are thermal,


where, by thermal, we refer to the explicitly T-dependent second
term in the propagators as defined in Appendix <ref>.  Since the
thermal part of all propagators puts the particle on mass-shell,
the last option yields a term depending on the product δ(k^2)
 δ((p'+k)^2 - m_f^2). Since the external fermion is on-shell,
p'^2 = m_f^2; hence this product of delta functions requires p'
· k = 0. But δ(k^2) implies k^0 = ±|k⃗|
or the angle between k and p' should satisfy cosθ_kp' = ±
p'^0/|p⃗'⃗| or |cosθ_kp'| > 1, which
is impossible, so there is no phase space available in the case that both
propagators are on-shell. Hence we only need to consider the contribution
of the first diagram when either one of the photon or the fermion propagator
is thermal. In particular, as discussed in Section <ref>, we only
need to consider the IR-finite G-photon contribution with the
g^μν term in the photon propagator being replaced by G^μν:

    -i
    g^μν   → -i  [ g^μν - b_k (p_f, p_i) k^μ k^ν]  , 
       ≡ -i  [ G_k^μν(p_f,p_i)] .

Since the virtual photon vertices in Diagram 1, Fig. <ref> are
on the final fermion leg and the scalar, we have p_f=p', p_i=p, with
b_k(p_f, p_i) defined in Eq. <ref>. We consider in turn, the
contribution from the thermal parts of the photon and fermion
propagators respectively to the cross section for the process shown
in Diagram 1 of Fig. <ref>. We present the purely t-channel
contribution in detail in order to highlight some technicalities, and
then go on to the other contributions.



  §.§.§ NLO thermal photon contribution to Diagram 1


The matrix element at NLO from the (t-channel) Diagram 1 in
Fig. <ref> when the G-photon propagator is thermal, is
given by,

    M^t_NLO ( Diagram  1, γ)     = 
    -∫d^4k/(2π)^4i e^2 λ ^2/2 k · p'(2πδ(k^2)n_B(| k^0|) ) Δ(l) Δ(l+k)
    
               [(u(p',m_f)γ_μ(k+p'+m_f) P_R  u(q',m_χ) )
    . 
                   .
    (k-2 p+2 q)_ν(v(q,m_χ) P_L   v(p,m_f) )]
     G_k^μν(p',p) , 
       ≡∫d^4k/(2π)^4(2πδ(k^2) n_B(| k^0 |) ) F_NLO^t,1γ(k) ,

where n_B is the Bose distribution function given in Eq. <ref>
and we have used the index ( Diagram  1,
γ) to indicate that this is the contribution from Diagram 1 of
Fig. <ref> when the photon is thermal.




At this point, we are in a position to justify the simplification of the
scalar propagator in the heavy-scalar approximation. The scalar propagator
for the t-channel diagrams contains either the inverse of (l^2 -
m_ϕ^2), or ((l+k)^2 - m_ϕ^2), with l ≡ q -p ≡
p'-q'. Since m_χ≪ m_ϕ, l^2 < m_ϕ^2. Since the photon
propagator is thermal, k^2=0 and l· k ∝| k^0
| is small since a large k^0 = |k⃗| will suppress
the contribution due to the presence of the n_B(| k^0 |)
term. Hence either of the scalar propagators can be approximated by just
the (m_ϕ^2)^-1 terms.

A similar reasoning will hold in the next Section <ref> when we
consider the fermion propagator to be thermal.  Then we can write
(l+k)^2 = ((p'+k)-q')^2 ≡ (t-q')^2 (see Eq. <ref>). Now
(p'+k)^2 ≡ t^2=m_f^2 since the fermion propagator is thermal
and t· q should be small, otherwise the contribution will be killed
by the n_F(| t^0|) term; see Eq. <ref>. A similar argument
holds for Diagrams 2 and 3 in Fig. <ref> when the anti-fermion
is thermal, where we can write (l+k)^2 = (q-p+k)^2 ≡ (q-t)^2. For
the u-channel matrix elements, the scalar momentum is l' = q'-p and
the same reasoning applies. Hence we will replace Δ(l) = i /(l^2 -
m_ϕ^2) → i/(-m_ϕ^2), etc., in this section. When the scalar
is heavy, but not much larger than the dark matter particle or the
scale of the interaction (√(s)), we will reinstate these terms
and recompute the cross section in Section <ref>.

The integration over the photon momentum k can be partly
completed using the delta function:

    ∫d^4k (2πδ(k^2)) F(k)     = 
    2π∫_-∞^∞dk^0 ∫_0^∞ K^2 dK
    ∫dΩ_k [δ((k^0)^2-K^2)] F(k^0, K,
    Ω_k) , 
        = 2 π∫dk^0 ∫dΩ_k
    ∫ K^2 dK [ δ(k^0 -K) + δ(k^0+K)
    ]/| 2k^0 | F(k^0, K, Ω_k) , 
        = π∫dΩ_k
    [∫_0^∞| k^0 |dk^0 F(k^0, k^0, Ω_k) +
    ∫_-∞^0 | k^0 |dk^0 F(k^0, -k^0, Ω_k) ] ,
    
       ≡π∫_0^∞ωdω[
    ∫dΩ_k F_+(ω, ω, Ω_k) + 
    ∫dΩ_k F_-(-ω, ω, Ω_k) ] ,

where K ≡|k⃗|. Hence there are two contributions
to M_NLO^t( Diagram 1,γ) (and to each
such matrix element), one where K → k^0 ≡ω, and the other
where K → -k^0 ≡ -ω, as can be seen from Eq. <ref>.
Note that the lower limit (ω→ 0) can be safely taken precisely
because the G photon contribution is guaranteed to be IR-finite.




The purely t-channel thermal photon contribution from Diagram 1 is
therefore given according to Eq. <ref> as

    σ_NLO^t( Diagram 1,γ)     = 1/32π s(2π)^4 |p⃗'⃗|/|q⃗'⃗|∫dcosθ[∑_spins(M_LO^t)^†M_NLO^t,1γ + h.c.
    ] , 
        = 1/32 s(2π)^4 |p⃗'⃗|/|q⃗'⃗|∫ωdω  n_B(ω) [
    ∫dcosθ∫dΩ_k
    [ F_+^t,1,γ + F_-^t,1,γ]
    ]  , 
       ≡1/32 s(2π)^4 |p⃗'⃗|/|q⃗'⃗|∫ωdω  n_B(ω)
    Int_ Diagram  1, γ^t .

On performing the angular integrations, the contribution of this term to the
NLO cross section is given by,

    Int^t_ Diagram  1,γ    = 
     64 π e^2 λ ^4/3 m_ϕ^6 P'(4 P'(3 H^4+P^2  P'^2) +
    logH-P'/H+P' 3H (H^2+P^2) (H^2+P'^2)
    ) .

Note that the result in Eq. <ref> is independent of ω;
in fact, the individual F_+ and F_- contributions to
Int_ Diagram  1, γ^t contain some
apparent divergent terms of the order of
1/ω. When the two terms are added, these cancel, leaving behind IR
finite terms that are integrable in ω. This is a reflection of
the fact that the G photon insertion was tailored precisely to yield
such an IR-finite result. The logarithmic terms arise
when p' or p are collinear with q' or q. These will drop out of
the final calculation[The collinear terms cancel between the G
and G contributions from the virtual and real corrections
respectively <cit.>.]

Similarly, the u-channel and crossed tu-channel NLO contributions
from Diagram 1 of Fig. <ref> when the photon propagator is thermal can be calculated from the corresponding u-channel NLO matrix
element which is given by

    M^u_NLO ( Diagram  1, γ)     = 
    - ∫d^4k/(2π)^4i e^2 λ^2/2 k · p'(2πδ(k^2)n_B(| k^0|) )
    Δ(l') Δ (l'+k) 
               [(u(p',m_f)γ_μ(k+p'+m_f)
    P_R   u(q,m_χ) ) (k-2 p+2 q')_ν. 
               . (v(q',m_χ)
     P_L   v(p,m_f) ) ] G_k^μν(p',p) .

Then the total NLO contribution from (the t-, u and cross tu
channels of) Diagram 1, Fig. <ref>, with thermal photon
propagator is given by,

    Int_ Diagram  1, γ^tt+uu-tu    =
    64 π e^2 λ ^4/3 m_ϕ^6 P'[24 H^4 P'-6 H^2 m_χ^2
    P' -6 m_χ^2 P' (m_f^2+3 P'^2)+8 P^2 P'^3
    . 
           . +logH-P'/H+P'(6 H^5+H^3 (6
    (P^2+P'^2)-3 m_χ^2)+
    3 H P'^2 (2 P^2-3 m_χ^2)) ] , 
        = 
    64 π e^2 λ ^4/3 m_ϕ^6 P'[4 P' (8 H^4-2  H^2 (4 m_χ^2+m_f^2)+
    5 m_χ^2 m_f^2) . 
       .      +3 logH-P'/H+P'(8 H^5-4 H^3 (2 m_χ^2+m_f^2)+5 H m_χ^2 m_f^2)
    ] , 
        non-coll⟶256 π e^2 λ ^4/3 m_ϕ^6(8 H^4-2  H^2 (4 m_χ^2+m_f^2)+
    5 m_χ^2 m_f^2) ,

where we have substituted for P'^2=H^2-m_f^2, P^2=H^2-m_χ^2
and dropped the collinear logarithmic terms in the last line. It can be
seen that Int_ Diagram 1,γ^tt+uu-tu given in
Eq. <ref> is independent of ω; hence the final
contribution from thermal photons to Diagram 1, Fig. <ref>
is

    σ_NLO^tt+uu-tu( Diagram  1, γ)     = 
     1/32 s(2π)^4 P'/P∫ωdω  n_B(ω)
    Int_ Diagram  1, γ^tt+uu-tu , 
        = 1/32 s(2π)^4 P'/P π^2 T^2/6×
    Int_ Diagram  1, γ^tt+uu-tu ,

where the T^2 temperature dependence of the cross section arises due
to the ω-independence of Int_ Diagram 1,γ^tt+uu-tu,
with 

    ∫_0^∞ωdω  n_B(ω)
    = π^2 T^2/6 .

We now discuss the contribution at NLO from Diagram 1,
Fig. <ref>, when the fermion propagator is thermal.



  §.§.§ NLO thermal fermion contribution to Diagram 1



As mentioned earlier, the thermal contributions from Diagram 1 only
arise when either the photon or fermion propagator is thermal,
that is, the explicitly T-dependent second term of the propagators in
Eqs. <ref> and <ref> contribute; the contribution
when both are thermal vanishes. We now consider the contribution
at NLO of the t-channel Diagram 1, Fig. <ref>, when the
photon propagator is non-thermal, but the fermion propagator is thermal. Recall that only the photon contributions give rise to IR
divergences (both at T=0 and finite temperature). This is because
of their differing thermal distributions; see Eqs. <ref> and
<ref>. In addition, we have shown in Section <ref>
that the K-photon contribution when the photon is inserted into a thermal fermion vanishes; see Eq. <ref>. Hence in
the G photon insertion, which arises from G^μν = (g^μν
- b_k k^μ k^ν), the second term, which is the K photon piece,
vanishes and we therefore need consider just the g^μν contribution
(or in other words, since the K contribution vanishes, the G term
contains the entire g^μν contribution). The corresponding
matrix element is then given by

    M^t_NLO ( Diagram  1, f)     =
    - ∫d^4k/(2π)^4i e^2 λ ^2   g^μν/k^2 (-2πδ((p'+k)^2-m_f^2)n_F(| (p'+k)^0|) )
    
               Δ(l) Δ(l+k)
    [(u(p',m_f)γ_μ(k+p'+m_f) P_R  u(q',m_χ) )
    . 
               .
    (k+2 p'-2 q')_ν(v(q,m_χ) P_L   v(p,m_f) )] ,

where n_F is the fermion number operator defined
in Eq. <ref>, the 1/k^2 term arises from the temperature-independent part of the photon propagator, and the
delta-function from the thermal fermion propagator. As before, the
index ( Diagram 1, f) refers to the contribution arising
from Diagram 1, Fig. <ref>, when the thermal part of the
fermion propagator contributes. The sign
difference between n_F and n_B is what dictates the IR finite nature
of the fermionic contributions. It is convenient to define p'+k ≡
t, and change the variable of integration to t, so that

    M^t_NLO ( Diagram  1, f)     =
    ∫d^4t/(2π)^4i e^2 λ ^2/2g^μν/m_f^2 - t · p' (2πδ(t^2-m_f^2)n_F(| t^0|) )
    
           Δ(l) Δ(l+k)
    [(u(p',m_f)γ_μ(t+m_f) P_R  u(q',m_χ) )
    . 
           .
    (t+p'-2 q')_ν(v(q,m_χ) P_L   v(p,m_f) )] ,

so that the same simplification can be done using the delta-function as
discussed in Eq. <ref>, to obtain 

    ∫d^4t [2πδ(t^2-m_f^2)]  F(t)     =
    2π∫d^4t  δ(t_0^2-|t⃗|^2 -m_f^2) F(t) 
    ≡ 2π∫d^4t  δ(t_0^2-ω_t^2) F(t) 
     , 
        = π∫_m_f^∞ K_t 
     dω_t [
    ∫dΩ_t F_+(ω_t, ω_t, Ω_t)
    
    + ∫dΩ_t F_-(-ω_t, ω_t, Ω_t) ],

where ω_t^2 = |t⃗|^2 + m_f^2 ≡ K_t^2 + m_f^2
and we have expressions for F_±(t) analogous to that in
Eq. <ref>. Then the thermal fermion contribution to the
purely t-channel Diagram 1 is given analogously to that for thermal photons in Eq. <ref> by

    σ_NLO^t( Diagram 1,f)     = 1/32π s(2π)^4 |p⃗'⃗|/|q⃗'⃗|∫dcosθ[∑_spins(M_LO^t)^†M_NLO^t,1 f + h.c.
    ] , 
        = 1/32 s(2π)^4 |p⃗'⃗|/|q⃗'⃗|∫ K_t dω_t   n_F(ω_t) [
    ∫dcosθ∫dΩ_t
    [ F_+^t,1,f + F_-^t,1,f]
    ]  , 
       ≡1/32 s(2π)^4 |p⃗'⃗|/|q⃗'⃗|∫ K_t dω_t   n_F(ω_t)
    Int_ Diagram  1, f^t .

Again, on completing the angular integration, we find the t-channel
thermal fermion contribution to Diagram 1, Fig. <ref>, to be

    Int_ Diagram  1, f^t    = 
    32 π e^2 λ ^4/3 K_t m_ϕ^6 P'[4 K_t P' (3 H^2 m_χ^2 +
    P^2 (4 H^2 - m_f^2) ) . 
            -(6 H^5 ω_t-3 H^4 m_f^2+ H^3 ω_t (4 P^2-3 m_f^2)
    -4 H^2 m_f^2 P^2-H ω_t P^2 (m_f^2-2 P'^2)
    . 
            - . 
    m_f^2 P^2 (P'^2 - m_f^2))
    logH ω_t+K_t P'-m_f^2/H ω_t-K_tP'-m_f^2
            +(6 H^5 ω_t+3 H^4 m_f^2+H^3 ω_t
    (4 P^2-3 m_f^2)
    +4 H^2 m_f^2 P^2-H ω_t P^2 (m_f^2-2 P'^2)
    . 
            + . . m_f^2 P^2 (P'^2-m_f^2))
    logH ω_t-K_t P'+m_f^2/H ω_t+K_t P'+m_f^2] .

While the non-logarithmic terms are (as in the case of thermal photons)
independent of ω_t, the logarithmic terms have a complicated
dependence on ω_t through K_t and it is not possible to
analytically integrate these terms. A simplification is achieved in
the limit that the fermion masses can be neglected in comparison
to ω_t, so that K_t →ω_t. Then the contribution
simplifies[Note that no approximations have been made in the
non-logarithmic term.] to

    Int_ Diagram  1, f^t   K_t →ω_t⟶64 π e^2 H λ ^4/3 m_ϕ^6 P'[0pt16pt 2 H P'(3m_χ^2 + 4 P^2) . 
               . - 2 (3 H^4 + 2 H^2 P^2 + P^2 P'^2)
    logH +P'/H-P'] ,

with

    log [H ω_t + K_t P' ± m_f^2]/log [H ω_t - K_t P' ± m_f^2]   K_t→ω_t⟶log [H + P']/log [H - P'] ,

where we have dropped terms of order O(m_f^2/(Hω_t)), so
that these logarithmic terms are of the same form as before, independent
of ω_t, yielding again a T^2 temperature dependence from this
contribution, using[Expanding K_t = √(ω_t^2 - m_f^2)
as K_t ≈ω_t - m_f^2/(2ω_t) + ⋯ is not a possible
choice since ∫(  dω_t /ω_t) n_F is not convergent;
it is of course possible to integrate this numerically if more precise
results are required.]

    ∫_m_f^∞ K_t dω_t   n_F(ω_t) ≈∫_0^∞ω_t dω_t   n_F(ω_t)
    = π^2 T^2/12 .

Similarly, the u-channel matrix element for thermal fermions is given
by:

    M^u_NLO ( Diagram  1, f)     =
    - ∫d^4t/(2π)^4i e^2 λ ^2
     g^μν/2(m_f^2 - t · p') (-2πδ(t^2-m_f^2)n_F(| t^0|) )
    
               Δ(l') Δ (l'+k)
    [(u(p',m_f)γ_μ(t+m_f) P_R  u(q,m_χ) ) .
    
               .
     (t+p'-2 q)_ν(v(q',m_χ) P_L   v(p,m_f) )] .

Using the LO matrix elements given in Eq. <ref>, the NLO matrix
element for the t-channel contribution given in Eqs. <ref>
and <ref> and that for the u-channel contribution in
Eq. <ref>, the total thermal fermion contribution from
Diagram 1, Fig. <ref>, that is, from the combined t-channel,
u-channel, and the crossed tu-channel, is found to be

    Int_ Diagram 1,f^tt+uu-tu    =
    32 π e^2 λ ^4/3 m_ϕ^6 K_t P'[4 K_t P'(3 m_χ^2 m_f^2 +
    2P^2(4 H^2 -m_f^2) ) . 
            +3 m_χ^2 logH ω_t + K_t P' - m_f^2/H ω_t - K_t P' - m_f^2( 3H^3ω_t - 2 H^2 m_f^2
    + H ω_t (P'^2 - 2 m_f^2) + m_f^4 - P'^2 m_f^2) 
            -3 m_χ^2 logH ω_t + K_t P' + m_f^2/H ω_t - K_t P' + m_f^2(3H^3ω_t + 2 H^2 m_f^2 + H
    ω_t (P'^2 - 2 m_f^2) - m_f^4 + P'^2 m_f^2) 
            -2 logH ω_t + K_t P' - m_f^2/H ω_t - K_t P' - m_f^2( 6H^5ω_t - 3 H^4 m_f^2 + H^3
    ω_t (4P^2 - 3 m_f^2) - 4H^2P^2 m_f^2 . 
                       . - H ω_t P^2(m_f^2 - 2
    P'^2) - P^2 m_f^2( P'^2- m_f^2)) 
            -2 logH ω_t + K_t P' + m_f^2/H ω_t - K_t P' + m_f^2(6H^5ω_t + 3 H^4 m_f^2 + H^3
    ω_t (4P^2 - 3 m_f^2) + 4H^2P^2 m_f^2 . 
                       . . - H ω_t P^2(m_f^2 - 2
    P'^2) + P^2 m_f^2( P'^2- m_f^2)) ] .

As in the case of the t-channel contribution alone, we can again use
the approximation given in Eq. <ref>, so that

    Int_ Diagram  1, f^tt+uu-tu   K_t→ω_t⟶128 π e^2 λ ^4/3 m_ϕ^6 [( 3 m_χ^2 m_f^2 + 2P^2(4 H^2 -m_f^2) )
    -3 H/2 P'logH+P'/H-P'{4 H^4+ . . 
                       . .
     2 H^2 (8 P^2-2 m_χ^2 - m_f^2) + m_f^2
    (3 m_χ^2 - 2 P^2 ) }0pt16pt] .

Again Int_NLO^1,tt+uu-tu,f is independent of ω_t; hence,
the ω_t integral (see Eq. <ref>) gives a
T^2 temperature dependence to the cross section (analogue of
Eq. <ref> for all channels) from the thermal fermions
in Diagram 1, Fig. <ref>, as well.



 §.§ Total thermal contribution to the NLO cross section


Now that we have demonstrated details of the calculation of the cross
section from Diagram 1, Fig. <ref>, we present the detailed
results from all Diagrams in Fig. <ref>. Diagram 3 is similar
to Diagram 1; here also, it is clear from an analysis similar to that
for Diagram 1 in Section <ref> that only one of the photon and
fermion propagators can be thermal at a time. Note that Diagrams 4
and 5 have no fermion propagators while Diagram 2 has three propagators
(excluding that of the scalar) and can in principle have contributions
from any one, any two, or all three propagators being thermal
(where we remind the reader that by thermal we mean the contribution
from the explicitly T-dependent second term in the expressions for
the photon and fermion propagators given in Eqs. <ref>
and <ref>). As argued for Diagrams 1 and 3, there is no
consistent available kinematic phase space when the photon propagator
and either one of the fermion or anti-fermion propagators are both
thermal. Hence there is no contribution from the case when all
three propagators in Diagram 2 are thermal (all three propagators
contribute on-shell) as well as when the photon propagator and either
the fermion or anti-fermion propagator is thermal. Left is the
case when the fermion and anti-fermion propagators are thermal,
but not the photon one. Here, we find that the delta function constraint
δ((p'+k)^2 - m_f^2) δ((p-k)^2 - m_f^2) is satisfied for the
single phase space point, k^0=0. Therefore, for all the diagrams, as
long as we discard the negligible contribution from scalar thermal
propagators, we only have thermal contributions when exactly one of the
propagators in the diagram is thermal and the others contribute
through their non-thermal or temperature-independent parts.

We now go on to present the contributions from the remaining Diagrams
2–5 in Fig. <ref> as well as their u-channel
counterparts. As can be observed from the diagrams in
Fig. <ref> itself, Diagrams 1, 3 and 5 have two scalar propagators,
Diagram 4 has three, and Diagram 2 as well as the LO matrix element have
one scalar propagator each. Hence it is expected that the contributions
of Diagram 4 will be suppressed compared to the other contributions.
As before, we have only one propagator whose thermal part contributes,
and we separately list the contributions from the terms when the photon
propagator is thermal and when the fermion/anti-fermion one is
thermal. We begin with the thermal photon contributions.



  §.§.§ Total thermal photon contributions to the cross section


The NLO matrix elements for the remaining Diagrams 2–5,
Fig. <ref>, when the thermal part of the photon propagator
is taken into account are analogous to the expression given in
Eq. <ref> and are given in Appendix  <ref>. The
corresponding contributions from Diagrams 2–5 to the cross section
(corresponding to Eqs. <ref>, <ref>,
<ref> from Diagram 1) from both the t- and u-channels
and crossed tu-channel (in the heavy scalar propagator limit)
are given by,

    Int_ Diagram 2,γ^tt+uu-tu    =
     64 π e^2 λ ^4/3 m_ϕ^4[
     (6 H^2-6 m_χ^2+2 P^2) + 
    1/HP'logH-P'/H+P'(3 H^2-P^2)  (H^2-P'^2) ] ,


    Int_ Diagram 3,γ^tt+uu-tu    =
     64 π e^2 λ ^4/3 m_ϕ^6 P'[2 P' (6 H^4+H^2 (8 P^2-3 m_χ^2)
     +(3 m_χ^2-2 P^2)
    (m_f^2-P'^2)) . 
            +H logH-P'/H+P'(12 H^4+H^2 (-9 m_χ^2-6 m_f^2+8 P^2) . 
           . . +m_χ^2
    (6 m_f^2-3 P'^2)-2 P^2 (m_f^2-2 P'^2)) ] ,


    Int_ Diagram 4,γ^tt+uu-tu    =
    512 π e^2 λ ^4/15 m_ϕ^8 P' (H^2-P'^2)[P' {60 H^8-30 H^6 (2 m_χ^2-m_f^2+2 P'^2)
    . . 
            +15 H^4 (m_χ^4-m_χ^2 (m_f^2-2 P'^2)
    +2 P^2 (3 m_f^2+2 P'^2))
            -H^2 [5 m_χ^2 (m_f^2 (P^2+3 P'^2)
    +2 P'^2 (P^2-3 P'^2) ) . 
           . +2 P^2 (m_f^2
    (15 P'^2-4 P^2)+30 P'^4)]
           . -P'^2 (15 m_χ^4 P'^2+5 m_χ^2 P^2
    (m_f^2-2 P'^2)+2 m_f^2 P^4)}
            + logH-P'/H+P' H m_f^2 (H^2-P'^2)
    {30 H^4-5 H^2 (3 m_χ^2-8 P^2) . 
           . . -15 m_χ^2 P'^2+4 P^4+10 P^2 P'^2}] ,


    Int_ Diagram 5,γ^tt+uu-tu    =
    512 π e^2 λ^4/3 m_ϕ^6(-6 H^4+3 H^2 m_χ^2+P'^2 (3 m_χ^2-2 P^2)) .

Simplifying, and substituting P'^2 → H^2 - m_f^2 and P^2 → H^2 -
m_χ^2, and dropping the collinear terms, the total contribution
from all terms where the virtual photon propagator is thermal, is given by

    Int_Total,γ T^tt+uu-tu    =
    {[Int_ Diagram 2,γ^tt+uu-tu] +
    [Int_ Diagram 1,γ^tt+uu-tu +
    Int_ Diagram 3,γ^tt+uu-tu +
    Int_ Diagram 5,γ^tt+uu-tu] . + 
       .         [Int_ Diagram 4,γ^tt+uu-tu]
    } , 
        = 128 π e^2 λ ^4/3{1/m_ϕ^4[4 P^2] - 2/m_ϕ^6 H^2 [H^2-m_χ^2-P^2 ]
    + 4/5 m_ϕ^8[90 H^6-120 H^4(m_χ^2-P^2). .
    
            + H^2(30 m_χ^4 + 5 m_χ^2(9 m_f^2 - 4 P^2)
    - 6 P^2(5 m_f^2 - P^2)) 
                   . . +m_f^2(-15 m_χ^4 +
     15 m_χ^2 P^2 + 2 P^4) ] 0mm8mm} , 
        = 512 π e^2 λ ^4/15{5/m_ϕ^4( H^2 - m_χ^2) + H^2/m_ϕ^8[
    216 H^4- 272 H^2 m_χ^2 + 56 m_χ^4 ] . 
       .              -m_f^2 /m_ϕ^8[ 28 H^4 - 86 H^2 m_χ^2 +
     28 m_χ^4 ] } ,

where we have combined contributions having the same number of
scalar propagators. A few points can be noted. From the first line,
we see that the total thermal photon contribution from Diagrams 1, 3,
and 5 vanishes, since P^2 = H^2 - m_χ^2. Hence terms
contributing as 1/m_ϕ^6 from thermal photon insertions vanish. Also,
from the last line in Eq. <ref>, we see that the contribution
from Diagram 2 (contributing as 1/m_ϕ^4) is proportional to H^2 -
m_χ^2 = P^2 and is small in the non-relativistic limit. So also
is the contribution from Diagram 4 (1/m_ϕ^8) term since the
term proportional to H^2 also vanishes when H = m_χ (exact
non-relativistic limit), leaving only the m_f^2 dependent term.
The familiar independence on ω is seen for the total thermal
photon contribution at NLO, so that the total cross section from these
terms again has a T^2 temperature dependence. We now calculate the
contribution from thermal fermions.



  §.§.§ Total thermal fermion contributions to the cross section


Only Diagrams 1–3 in Fig. <ref> contribute to the NLO cross
section when the thermal part of the fermion propagator is considered. The
corresponding matrix elements for Diagrams 2 and 3 are also given in
Appendix <ref>. The contributions to the NLO cross section from
each of these diagrams (in the heavy scalar limit, where we
have substituted for P' and dropped the logarithmic terms, which have
the same form as for Diagram 1), are given by

    Int_ Diagram 2, f^tt+uu-tu    =
    -32 π e^2 λ^4/3 m_ϕ^4 (H^2 - ω_t^2)(12 H^2(H^2 - m_χ^2) +
    4 H^2 P^2 + 2 ω_t^2 (3 (H^2-m_χ^2) + P^2) ) , 
    
    Int_ Diagram 2, f^tt+uu-tu    =
    Int_ Diagram 2, f^tt+uu-tu , 
    
    Int_ Diagram 3,f^tt+uu-tu    =
    128 π e^2 λ^4/3 m_ϕ^6( 3 m_χ^2m_f^2 + 2 P^2(4H^2- m_f^2) ) .

Notice that the contribution of the thermal fermion from Diagram
3 is the same as that from Diagram 1, as can be seen[Note the
absence of odd powers of ω_t due to the symmetry explicit in
Eq. <ref>.] from Eq. <ref>. Furthermore, the contribution of
the thermal fermion from Diagram 2, Fig. <ref> has non-trivial
ω_t dependences both in the numerator and denominator. Recognising
that large values of ω_t will lead to vanishing of the
corresponding distribution function, n_F(|ω_t|), we
expand the denominator for ω_t < H (recall that H ≥ m_χ and
β m_χ = x ∼ 20 at freeze-out) to obtain the total thermal fermion contribution from Diagram 2 to be

    Int_ Diagram 2, f+f^tt+uu-tu    =
    -64 π e^2 λ^4/3 m_ϕ^4 H^2(12 H^2(H^2 - m_χ^2) +
    4 H^2 P^2 + 6 ω_t^2 (3 (H^2-m_χ^2) + P^2) ) .

The presence of the ω_t^2 terms in the numerator will lead to
T^4 temperature dependence in the cross section[Since
H^2 - m_χ^2 = p^2, both the ω_t^0 and the ω_t^2 terms
are small in the non-relativistic limit.]. Then, the total contribution
from all terms where the virtual fermion (or anti-fermion) propagator
is thermal, is given by

    Int_Total,(f+f) T^tt+uu-tu    =
    [Int_ Diagram 1, f^tt+uu-tu +
    Int_ Diagram 3, f^tt+uu-tu] +
    [Int_ Diagram 2, f^tt+uu-tu +
    Int_ Diagram 2, f^tt+uu-tu]
     , 
        = 64 π e^2 λ^4/3 m_ϕ^6{[0pt14pt 12 m_χ^2m_f^2 + 8 P^2(4H^2- m_f^2) ] -m_ϕ^2
     [ 0pt14pt 12 (H^2 - m_χ^2) + 4 P^2 .. 
               . . + 6 ω_t^2/H^2(3 (H^2-m_χ^2) + P^2 ) ] } .

The total invariant NLO cross section is then given by the sum of the
thermal photon and thermal fermion contributions, 

    σ_NLO(s)     = 
    1/32 s(2π)^4 P'/P[∫ω dω  n_B Int_Total,γ T^tt+uu-tu + ∫ K_t  dω_t   n_F
    Int_Total,(f+f) T^tt+uu-tu] ,
    
        = λ^4/12 π s P'/P [8 α
    T^2/15] ×
     {0mm6mm1/m_ϕ^4[ 7 π^2 T^2/4 H^2 (H^2 - m_χ^2) ] +
    1/m_ϕ^6[ 0pt12pt 5 m_χ^2m_f^2 + 2 P^2(4H^2-
     m_f^2) ] . 
           .  + 1/m_ϕ^8[ 0pt14pt 8 H^2
    (27 H^4- 34 H^2 m_χ^2 + 7 m_χ^4) - 2 m_f^2 (14 H^4 - 43 H^2 m_χ^2 +
     14 m_χ^4) ] } ;

here s = 4H^2 and the terms have been ordered in increasing powers of
m_ϕ^2 in the denominator. We have replaced K_t by ω_t in
the integrand of the second term in order to obtain an analytical result
and furthermore, presented the results, taking the lower limit of the
integration to be zero. Notice that the leading term (1/m_ϕ^4) is
of order O(T^4) and is further suppressed by (H^2 - m_χ^2)=
P^2 in the non-relativistic limit. The terms suppressed by an additional
propagator factor (∼ 1/m_ϕ^6) contribute at T^2 order and are
either proportional to m_f^2 (∼ 5 m_χ^2 m_f^2) or to P^2,
which is small is the non-relativistic limit. However, these terms,
of order O(H^2 P^2 T^2/m_ϕ^6) or O(m_χ^2 m_f^2
T^2/m_ϕ^6), can have significant contribution in the relativistic
regime when P is large.



  
The NLO cross section in the non-relativistic limit:
Freeze-out occurs around m_χ/T ∼ 20 so that the dark matter
particles can be considered to be non-relativistic at this point,
with P→ m_χ v, and H^2 → m_χ^2(1 + v^2), so that
Eqs. <ref> and <ref> reduce to

    Int_Total,γ T^tt+uu-tu   v small⟶512 π e^2 λ ^4/15 m_ϕ^8{
     30 m_χ^4 m_f^2 (1+ v^2) + 5 m_ϕ^4 m_χ^2 v^2 } ,
     
    
    Int_Total,(f+f) T^tt+uu-tu   v small⟶64 π e^2 λ^4/3 m_ϕ^6{[ 6 m_χ^2 m_f^2 + 4 m_χ^2 v^2 (4 m_χ^2 - m_f^2) ]
    + 8 m_ϕ^2 m_χ^2 v^2 (3 ω_t^2/H^2 -
    2 ) } .

It can be seen that the NLO cross section (obtained from
Eqs. <ref> and <ref>) is proportional to m_f^2
as v→ 0, just as the LO cross section; see Eq. <ref>. This can
be understood from helicity conservation: just as in the LO case, the NLO
diagrams (with an additional virtual photon) that we have computed
are planar 2→ 2 processes. The Majorana coupling then forces the
final states into the “wrong” chirality so that the cross section is
proportional to the fermion mass squared. At early times in the evolution
of the Universe, the dark matter particles are relativistic, √(s)
can be large, and there is less suppression in the annihilation into the
lighter fermions. This also holds for freeze-in scenarios where m_χ/T
∼O(1). Now that we have understood the structure of the NLO
contributions in the heavy-scalar limit, we will go on to include the
entire scalar propagator in the calculation of the NLO cross section.



§ THE DARK MATTER ANNIHILATION CROSS SECTION AT
NLO INCLUDING THE “DYNAMICAL” SCALAR PROPAGATOR



We repeat the calculation of the previous section without making the
heavy scalar assumption, i.e., retaining the “dynamical”
scalar propagator, iD_ϕ = i/(l^2 - m_ϕ^2), where l =
q -p is the momentum of the scalar. As we had already argued in Section
<ref>, l · k, k^2 ≪ m_ϕ^2; hence, we expand the relevant
propagator terms, for instance,

    1/(l+k)^2 - m_ϕ^2    =
    1/(l^2 - m_ϕ^2) + (2 l· k + k^2)
       ≈1/(l^2 - m_ϕ^2)[1 - (2 l· k + k^2)/(l^2 - m_ϕ^2)] .

Similar expansions can be obtained for the u-channel scalar propagators
as well. This enables us to get analytic expressions for the NLO cross
section. However, the results are rather cumbersome and opaque; they are
available as Mathematica <cit.> Notebooks in Ref. <cit.>. We simply make
a few remarks on the calculation and specifically on the non-relativistic
limit which is of interest at freeze-out.

The main difference on including the momentum dependence of the scalar
is the increase in complexity of the cosθ angular integration
due to the presence of the term l^2 - m_ϕ^2 ≡ m_χ^2 + m_f^2 -
m_ϕ^2 - 2 q · p ≡ m_Φ^2 - 2(H^2 - P P' cosθ) in
the denominator. It will be seen that this gives rise to logarithms of
the type log(2 H^2-m_Φ^2 ± 2 P P') which we have
already encountered when calculating the LO cross section.

Once again, the leading O(T^2) terms are ω-
(ω_t)-independent for the thermal photon (fermion) cases, while
the ω^2 (ω_t^2) terms contribute to the O(T^4)
terms. We have,

    σ_NLO(s)     =   1/32 s (2π)^4P'/P[ ∫ dωω  n_B 
      Int_Total,γ T^tt+uu-tu +
    ∫ dω_t ω_t   n_F 
      Int_Total,(f+f) T^tt+uu-tu] ,
    
       v small⟶λ^4/4 π sP'/Pπα T^2/6[Int_NLO^a + v^2   Int_NLO^b ] ,

where the expression in the second line refers to the so-called a and b
terms that contribute in the non-relativistic limit.

It is interesting to note that the pattern is the same as with the heavy
scalar computation in the previous section. The leading behaviour from
naive counting of the scalar propagator is expected to have the form
(from Diagram 2)

    σ_NLO(s) ∝m_χ^2 T^2/(2 H^2 - m_Φ^2 + 2 P P') (2 H^2 - m_Φ^2 - 2 P P') ,

where the denominator (proportional to m_ϕ^4 in the heavy scalar
limit) is the analogous term to the 1/m_ϕ^4 contribution in the
heavy-scalar limit. However, as in the previous section with heavy-scalar
assumption, we find this term vanishes and the leading contribution
is again

    σ_NLO(s)    ∝m_ϕ^2 M^2 m_f^2 T^2/(4 H^4 - 4 m_Φ^2 H^2 + m_Φ^4 - 4 P^2 P'^2)^2 , 
       ∝M^2 m_f^2 T^2/m_ϕ^6 ,

where M^2 is any combination of H ≡√(s)/2, m_χ or
m_f, of dimension 2.

As mentioned earlier, the results for the general relativistic case are
available as sets of Mathematica Notebooks at the web-page
Ref. <cit.>. Here we will discuss some aspects of the results in
the non-relativistic limit, which are of relevance near freeze-out, when
m_χ/T ∼ 20.



  
NLO cross section in the non-relativistic limit:

In the limit when p → m_χ v, and the expressions are expanded
for small velocities, the cross section can be expressed as σ  v = a + b v^2. The contributions to the “a” term from all the
diagrams (see Eq. <ref> for the definition) 
are listed in Table <ref>. Both O(T^2)
and O(T^4) terms contribute. (The O(v^2) terms can
be found from the expressions given online in Ref. <cit.>.) In the
non-relativistic limit, the denominator in the full expressions reduces
to the form visible in Table <ref>:

    1/(4 H^4 - 4 H^2 m_Φ^2 + m_Φ^4 - 4 P^2 P'^2)^n   ≈1/(m_χ^2 - m_f^2 + m_ϕ^2)^n[ 1 - 4 n P^2
    m_ϕ^2/(m_χ^2 - m_f^2 + m_ϕ^2)] , 
    
       v→ 0⟶    1/(m_χ^2 - m_f^2 + m_ϕ^2)^n≡1/D^n ,

as was encountered in the LO expression for the cross section in
Eq. <ref>. It can be seen from Table <ref> that,
in the non-relativistic limit, the “a” terms in the cross section
are proportional to the square of the fermion mass, both at order
O(T^2) and O(T^4), just as in the LO case. Again,
as with the NLO calculation in the heavy-scalar limit, the potential
leading contribution at order 1/m_ϕ^4 is suppressed by a factor of
P^2 and contributes[Expressions for the “b” term in the cross
section are given online in Ref. <cit.>; note that these are not
necessarily proportional to m_f^2.] at the order O(P^2
T^4 m_ϕ^6/(m_χ^2 D^5)) ∼O(P^2 T^4 /(m_χ^2
m_ϕ^4)), as before. Hence, the leading contribution is again at order
O(1/m_ϕ^6), as mentioned earlier. Note that only Diagrams
2 and 4 have O(T^4) contributions; that from Diagram 4 is
highly suppressed due to the 1/D^5 ∼ 1/m_ϕ^10 denominator
factor. In general, it can be seen that the leading contributions are
of order O(m_f^2 m_χ^2 T^2/m_ϕ^6) in the
non-relativistic limit.



Considering only the corresponding “a” (v^0) terms in the
non-relativistic limit, the relative size of the NLO contribution for each flavour of fermion pair is given by,

    σ_NLO^a/σ_LO^a    = πα T^2/6
    m_ϕ^2m_f^2 (22 m_χ^2 + 3 m_f^2)/m_f^2 m_χ^2  ,
    
       ≈11 πα/3T^2/m_ϕ^2 ,

where we have used D ∼ m_ϕ^2 and dropped terms in
σ^NLO_a having higher powers of D=m_ϕ^2 in the denominator.



§ DISCUSSIONS AND CONCLUSION


We have computed the next-to-leading order thermal contributions to the
cross section for dark matter annihilation into fermion pairs. We have
used a Lagrangian which couples the dark matter Majorana fermions
to standard model fermions via charged scalar doublets. We use the real
time formulation of thermal field theory to compute the higher order
contributions to the cross section for the annihilation of dark
matter particles[Since the dark matter particles are taken
to be Majorana fermions, the annihilation occurs via the t-, u-
and cross tu-channel processes; Dirac dark matter particles will have
the same interaction where the tu cross terms vanish.] via χχ→ f f in a heat bath of fermions, scalars and photons at
temperature T. Note that the regime of interest is during the evolution
of the Universe when electroweak symmetry breaking is complete so that
we consider higher order corrections only from photons, since the W
and Z are by now massive. We assume in addition that the scalars are
heavy, m_ϕ≳ m_χ.

We have used the generalised Grammer and Yennie
<cit.> technique where
the cancellation of infra-red soft divergences is straightforward. This
is achieved by separating the photon propagator into the sum of K
and G parts. This can be realised in the thermal case as well,
where the propagators are now sums of temperature-independent and
temperature-dependent 2 × 2 matrices; we label the latter as
the thermal part of the propagator since it carries an explicit
temperature dependence. It was shown earlier <cit.> that the infrared (IR) divergences (which are much more
acute than at zero temperature, being linear rather than logarithmically
divergent in the soft limit) are completely contained in the K photon
contribution and the G photon term is finite. The infrared divergences
cancel against similar terms from real soft photon contributions; note
that both absorption from, and emission into the heat bath is required
to be included for these cancellations to occur.

While the use of the Grammer and Yennie technique (and its modification for
thermal field theories) is well-known for demonstration of the
cancellations of IR divergences in various contexts, it has  been rarely
used, so far as the authors are aware, to actually compute the IR-finite
remainder. In the present work, we have computed the NLO finite
G photon contributions to the dark matter annihilation cross section
using this technique and found its application to be straightforward. A
further simplification occurred since kinematically it turns out that the
total contribution to the NLO cross section is actually a sum of terms
where any one of the (photon, fermion, anti-fermion) propagators can
be thermal at a time. Hence the Grammer and Yennie approach allows
a great simplification of the problem and consequent ease of
computation.

The calculation was first done by replacing the full scalar propagator
1/(l^2 - m_ϕ^2) by just the mass term, 1/(-m_ϕ^2), valid when
the scalar is very heavy. This enabled us to have simpler expressions so
that the pattern of dependence of the NLO cross section on both the
temperature T as well as the energy scales, H = √(s)/2 and m_f,
was visible. The five diagrams that contribute at NLO are shown in
Fig. <ref> for the t-channel processes with analogous diagrams
for the u-channel ones.

The potential leading contribution[This is the order
at which the LO cross section contributes; see Eqs. <ref>,
<ref> and <ref>.] at order 1/m_ϕ^4 in
the scalar mass (from Diagram 2) is of order O(p^2T^4/(H^2
m_ϕ^4)), which is small in the non-relativistic limit when p^2 is
small. Hence the leading contribution is at order 1/m_ϕ^6. It was
found that the O(T^2) contribution to the NLO cross section
at this order is proportional to either (m_f^2 M^2 T^2/m_ϕ^6) or
to (p^2 M^2 T^2/m_ϕ^6) where M^2 ∼ H^2 or m_χ^2. Hence
in the non-relativistic regime (at around m_χ/T ∼ 20) where
freeze-out occurs, the momentum-dependent term can be ignored and the
NLO cross section is proportional to the square of the fermion mass;
this dependence is also seen in the LO cross section and occurs due to
the Majorana nature of the dark matter particles and consequent helicity
suppression. It is likely that the NLO terms retain this feature due to
the 2 → 2 planar nature of the contributing diagrams.

The calculation was repeated, retaining the full dynamical scalar
propagator. The essential features of the calculation remain the same.
Since the expressions are very long, we present here only the leading
“a” term of the cross section in the non-relativistic limit
(contribution to σ v when v→ 0) in Table <ref>; the
entire solution in the relativistic case as well as their non-relativistic
approximations are listed in the Mathematica Notebooks available
online <cit.>. We find that both the O(T^2) and the
O(T^4) contributions to the “a” term are proportional
to m_f^2, as can be seen from Table <ref>. As with the
heavy-scalar approximation, the potential leading contribution at order
1/m_ϕ^4 is suppressed by a factor of p^2 and contributes at the
order O(p^2 T^4 m_ϕ^6/(m_χ^2 D^5)) which reduces, in
the limit m_ϕ > m_χ, to the order O(p^2 T^4 /(m_χ^2
m_ϕ^4)), as before. Hence, the leading contribution is again at order
1/m_ϕ^6, leading to a relative NLO thermal contribution which is
∼ 10  (α T^2/m_ϕ^2) times the LO cross section.

In summary, the NLO thermal correction terms have a T^2 dependence,
but are suppressed by additional powers of the (heavy scalar) propagator
as well as by the square of the fermion masses. This was discussed
in Ref. <cit.> where the NLO thermal cross section was
initially computed to be of order O(T^2), and later it was
shown that the T^2 contributions vanish and the leading contribution
is in fact of order O(T^4). This result was strengthened
in Ref. <cit.> where an operator product expansion
(OPE) approach was used to show the same result (for Dirac-type
dark matter particles); here, the scalar propagator was
again reduced to just the mass term (1/(-m_ϕ^2)), and it was
argued that O(m_f^2 T^2/m_χ^4) are small compared to
O(T^4/m_χ^4) and were thus ignored.

In this work, we have retained the fermion masses and found that, in
the non-relativistic limit, the leading terms as v → 0 at both order
O(T^2) and O(T^4) are suppressed by m_f^2, just as
the LO cross section. The general results when T may be large, so x =
m_χ/T ∼O(1), are given without making the non-relativistic
approximation, in Mathematica notebooks <cit.>. These corrections
may be significant in a freeze-in scenario in the early Universe when T
is larger for light dark matter candidates where m_χ/T ≳ 1.

Thermal corrections to such cross sections can become important in the
early Universe where the thermally averaged cross section, ⟨σ v ⟩, appears as the collision term in the Boltzmann
equation for the evolution of the dark matter phase space densities,
or equivalently, their yields. We have

    ⟨σ· v_Mø l⟩≡∫σ· v_Mø lexp[-E_1/T] exp[-E_2/T]  d^3 p_1
     d^3 p_2/∫exp[-E_1/T] exp[-E_2/T]  d^3 p_1  d^3 p_2 .

Here p_1, p_2 are the momenta of the annihilating dark matter
particles, and v_Mø l = √((p_1 · p_2)^2 - m_χ^2).
The Boltzmann approximation for the dark matter number densities has been
used. For cross sections of the form σ v = a + b v^2, we have

    ⟨σ v ⟩ = a + 3/2 b T .

Hence, the thermal corrections that we have computed at NLO will yield a
correction,

    ⟨σ_NLO v ⟩    = ⟨ (a_1 T^2 + a_2 T^4) + v^2 (b_1 T^2
    + b_2 T^4) ⟩ , 
        = (a_1 T^2 + a_2 T^4) + 3/2 (b_1 T^3 + b_2 T^5) .

The corrections from the NLO thermal contribution to “b” (the terms
proportional to v^2 and higher powers) can therefore be neglected.
These thermal corrections to the cross section σ(s) can alter the
yield equation and hence the final relic density, especially in freeze-in
scenarios. This is acquiring importance <cit.> in
view of the increasingly precise measurements of the relic density
<cit.>. It will be interesting to
pursue this further, within the frame-work of the Grammer and Yennie
formulation where the calculation is simpler; this is beyond the scope
of the present work.







§ FEYNMAN RULES IN THERMAL FIELD THEORY



The scalar propagator is given by

    i S_ scalar^t_a, t_b (p,m) = (
    	[   Δ(p)      0;      0 Δ^*(p) ]) +
    	2 πδ(p^2-m^2) n_ B(| p^0 |) 
    	([               1 e^| p^0 | /(2T); e^| p^0 | /(2T)               1 ]) ,

where Δ(p) = i/(p^2-m^2+iϵ), and t_a, t_b (=1,2) refer
to the field's thermal type. We shall refer to the first term as the
temperature-independent part and the second as the thermal
part in the main text since it carries explicit temperature dependence;
note that the latter contributes only on mass-shell. Only the thermal parts can convert type-1 to type-2 fields, and vice versa;
type-1 fields are the physical fields.

In the Feynman gauge, the photon propagator corresponding to
a momentum k is given by

    i  D^t_a,t_b_μν (k) =
    -g_μν i D^t_a,t_b (k) =
    - g_μν  i S_ scalar^t_a, t_b (k,0) ,

and the fermion propagator at zero chemical potential is given by

    i S_ fermion^t_a, t_b (p,m)     = (
    	[   S   0;   0 S^* ])
    	- 2 π S' δ(p^2-m^2) n_ F(| p^0 |) 
    	([                       1  ϵ(p_0) e^| p^0 | /(2T); -ϵ(p_0) e^| p^0 | /(2T)                       1 ]) , 
       ≡ (p+m) ([  F_p^-1  G_p^-1; -G_p^-1 F_p^*-1 ]) , 
       ≡ (p+m) S^t_a,t_b(p,m) ,

where S = i/(p -m+iϵ), and S' = (p +m);
hence the entire fermion propagator is proportional to (p+m),
just as at T=0. It can be seen that all propagators are a sum of temperature-independent and thermal parts. The fermionic
number operator,

    n_ F(|p^0|) ≡1/exp{| p^0 |/T} + 1 p^0 → 0⟶ 1/2 ,

is well-defined in the soft limit; however, the bosonic number operator
contributes an additional power of k^0 in the denominator to the
photon propagator in the soft limit:

    n_ B(| k^0 |) ≡1/exp{| k^0 |/T} - 1 k^0 → 0⟶ T/| k^0 | ,

so that the leading IR divergence in the finite temperature part is
linear rather than logarithmic as was the case at zero temperature. (The
cancellation of IR divergence thus involves cancellation of the leading
linear divergence as well as the logarithmic sub-divergence.)

The fermion–photon vertex factor is given by
(-ieγ_μ)(-1)^t_μ+1, where t_μ=1,2 for the
type-1 and type-2 vertices.
The scalar–photon vertex factor is [-ie(p_μ
+ p'_μ)](-1)^t_μ+1 where p_μ (p'_μ) is the 4-momentum of
the scalar entering (leaving) the vertex, while the 2-scalar–2-photon
seagull vertex factor (see Fig. <ref>) is
[+2ie^2 g_μν](-1)^t_μ+1 (the factor `2' is dropped for a
tadpole vertex). At a given vertex, all fields are of the same
type.



The bino-scalar-fermion vertex factor is i λ P_L; for details
on Feynman rules for Majorana particles at zero temperature, see
Ref. <cit.>. An overall negative sign applies as usual to
the type-2 bino vertex; again all fields at a vertex are of the same type.



 §.§ Some identities at finite temperature



Various identities useful for fermions are given in Ref. <cit.>
and are reproduced here for completeness. Note that
Eqs. <ref>, <ref> and <ref> lead to

    (p-m)   iS_p^t_a,t_b    = i (-1)^t_a+1δ_t_a,t_b ,

where we have used the compressed notation,
iS^t_a,t_b(p,m) ≡ i S^t_a,t_b_p.

Consider the insertion of the μ vertex of the additional virtual
photon with momentum k between vertices μ_q+1
and μ_q on the p' fermion leg; see Fig. <ref>.
The momentum of the photon at the
vertex μ_q is l_q, with
Lorentz index μ_q, and thermal
type-index t_q. Hence the momentum of the fermion leg to the left
of the vertex μ_q is
p'+∑_i=1^ql_i which we denote as
p'+∑_q. Using Eq. <ref>, we have,

    S^t_q,t_μ_p'+∑_q k 
    	S^t_μ,t_q+1_p'+∑_q+k    =
    	(-1)^t_μ+1[S^t_q,t_q+1_p'+∑_qδ_t_μ,t_q+1 - S^t_q,t_q+1_p'+∑_q+kδ_t_μ,t_q] .

If the photon vertex is inserted to the right of the vertex labelled
`1' on the fermion leg with momentum p', we have,

    u(p') k S^t_μ,t_1_p'+k    =
    u(p') (-1)^t_μ+1δ_t_μ,t_1 ,

since p'u(p') = m u(p'). Similar relations hold for the
insertion of of the virtual K photon at a vertex ν on the (anti-)fermion
p leg since p  u(p) = m  u(p) as well.



§ MATRIX ELEMENTS FOR ΧΧ→ F F AT NLO



As described in the text, the matrix elements are sums of the
contributions when any one of the photon, fermion or anti-fermion
propagators is thermal, that is, contributes via its explicitly
temperature-dependent part. The t- and u- channel matrix elements
for Diagram 1 when either the photon or fermion propagator is thermal is given in the text, in Eqs. <ref>, <ref>,
<ref>, and <ref>. The remaining matrix elements are
listed below. In the equations that follow, the indices
(Diagram i, γ),
(Diagram i, f) and
(Diagram i, f) refer to the contribution from
Diagram i (i = 1–5), when the photon, fermion, or anti-fermion
respectively contribute via the explicitly thermal parts of
their propagators.

The thermal photon contributions from the last four diagrams
in Fig. <ref> are given by

    M^t_NLO ( Diagram  2, γ)
        = 
    ∫d^4k/(2π)^4i e^2 λ ^2/4 k · p' k· p(2πδ(k^2) n_B(| k^0 |) ) (i Δ(l+k)) 
           [(u(p',m_f)γ_μ (k+p'+m_f) P_R  u(q',m_χ) )
    . 
           . (v(q,m_χ) P_L
    (k-p+m_f) γ_ν  v(p,m_f) )
    ] G_k^μν(p',p) , 
    M^t_NLO ( Diagram  3, γ)     = 
    ∫d^4k/(2π)^4i e^2 λ ^2/2 k· p(2πδ(k^2) n_B(| k^0 |) ) Δ(l) Δ(l+k) 
           [(u(p',m_f) P_R  u(q',m_χ) )
     (k - 2p + 2 q)_μ. 
       .     (v(q,m_χ) P_L
    (k-p+m_f) γ_ν  v(p,m_f) )
    ] G_k^μν(p,p) , 
    M^t_NLO ( Diagram  4, γ)     =
    ∫d^4k/(2π)^4 (i e^2 λ ^2)
    (2πδ(k^2) n_B(| k^0 |) ) Δ(l) i Δ(l+k) Δ(l)
    
           [(u(p',m_f) P_R  u(q',m_χ) )
     (k - 2p + 2 q)_μ  (k - 2p + 2 q)_ν. 
           . (v(q,m_χ) P_L  v(p,m_f) )
    ] G_k^μν(p,p) , 
    M^t_NLO ( Diagram  5, γ)     = 
    ∫d^4k/(2π)^4
    (i e^2 λ ^2)
    (2πδ(k^2) n_B(| k^0 |) ) Δ(l) Δ(l) 
           [(u(p',m_f) P_R  u(q',m_χ) )
     (v(q,m_χ) P_L  v(p,m_f) )
    ] G_k^μν(p,p) ,

for the t-channel, and appropriately crossed ones for the u-channel,
for example, as shown in Eq. <ref> for the contribution from
the first diagram. The matrix elements for the case when the fermion or
anti-fermion is thermal arise only from Diagrams 1, 2, and 3. They
are given by

    M^t_NLO ( Diagram  2, f)     =
    ∫d^4t/(2π)^4i e^2 λ ^2
     g^μν/4 (m_f^2-p' · t) (p· p'+m_f^2 - t ·
    (p+p'))
           (-2πδ(t^2-m_f^2)n_F(| t^0|) )
    i Δ(l+k) [(u(p',m_f)γ_μ(t+m_f) . . 
       . .      P_R  u(q',m_χ) )
    (v(q,m_χ) P_L  (t-p' - p) γ_ν
    v(p,m_f) ) ] ,


    M^t_NLO ( Diagram  2, f)     =
    ∫d^4t/(2π)^4i e^2 λ ^2
     g^μν/4(m_f^2-p · t) (p· p'+m_f^2 - t ·
    (p+p'))
           (-2πδ(t^2-m_f^2)n_F(| t^0|) )
     i Δ(l+k) [(u(p',m_f)γ_μ(-t+p' + p+m_f)
    . . 
       . .      P_R  u(q',m_χ) )
     (v(q,m_χ) P_L  (-t + m_f ) γ_ν
    v(p,m_f) ) ] , 
    M^t_NLO ( Diagram  3, f)     =
    ∫d^4t/(2π)^4i e^2 λ ^2
     g^μν/2(m_f^2 - t · p)(-2πδ(t^2-m_f^2)n_F(| t^0|) ) Δ(l)
     Δ(l+k) 
           [(u(p',m_f)
    P_R  u(q',m_χ) )
     (-t-p+2 q)_μ. 
           . (v(q,m_χ) P_L  (-t+m_f ) v(p,m_f) ) ] .

Here f corresponds to the fermion being thermal in Diagram 2 and
f corresponds to the anti-fermion being thermal in Diagrams
2 and 3 respectively. We use the same replacement technique here as
explained in Eqs. <ref> and <ref> for the thermal
fermion contribution to Diagram 1. For the matrix element of Diagram 2,
we have substituted t = p'+k as for Diagram 1; and we have substituted
t=p-k when the anti-fermion is thermal. Analogous expressions hold
for the u-channel matrix elements.



