Dark Matter Sommerfeldenhanced annihilation
and Boundstate decay at finite temperature
Abstract
Traditional computations of the dark matter (DM) relic abundance, for models where attractive selfinteractions are mediated by light forcecarriers and bound states exist, rely on the solution of a coupled system of classical onshell Boltzmann equations. This idealized description misses important thermal effects caused by the tight coupling among forcecarriers and other charged relativistic species potentially present during the chemical decoupling process. We develop for the first time a comprehensive abinitio derivation for the description of DM longrange interactions in the presence of a hot and dense plasma background directly from nonequilibrium quantum field theory. Our results clarify a few conceptional aspects of the derivation and show that under certain conditions the finite temperature effects can lead to sizable modifications of DM Sommerfeldenhanced annihilation and boundstate decay rates. In particular, the scattering and bound states get strongly mixed in the thermal plasma environment, representing a characteristic difference from a pure vacuum theory computation. The main result of this work is a novel differential equation for the DM number density, written down in a form which is manifestly independent under the choice of what one would interpret as a bound or a scattering state at finite temperature. The collision term, unifying the description of annihilation and boundstate decay, turns out to have in general a nonquadratic dependence on the DM number density. This generalizes the form of the conventional LeeWeinberg equation which is typically adopted to describe the freezeout process. We prove that our number density equation is consistent with previous literature results under certain limits. In the limit of vanishing finite temperature corrections our central equation is fully compatible with the classical onshell Boltzmann equation treatment. So far, finite temperature corrected annihilation rates for longrange force systems have been estimated from a method relying on linear response theory. We prove consistency between the latter and our method in the linear regime close to chemical equilibrium.
Contents
 I Introduction
 II Realtime formalism prerequisites
 III Equations of motion in Realtime formalism for nonrelativistic particles in a relativistic plasma environment
 IV Twotime BetheSalpeter equations
 V Twoparticle spectrum at finite temperature

VI DM number density equation in grand canonical ensemble
 VI.1 Recovering the LeeWeinberg equation
 VI.2 Spectral function, Sommerfeld enhancement factor and decay width for vanishing thermal corrections
 VI.3 DM number density equation for vanishing thermal corrections
 VI.4 Chemical potential in narrow thermal width approximation
 VI.5 Comparison to linearresponsetheory method
 VII Numerical results for twoparticle spectral function at finite temperature
 VIII Discussion
 IX Summary and conclusion
 A Semigroup property of free correlators
 B Annihilation term
 C Hard thermal loop approximation
 D Vacuum limit of twoparticle spectral function
 E Number density and chemical potential in grand canonical ensemble
I Introduction
The cosmological standard model successfully describes the evolution of the largescale structure of our Universe. It requires the existence of a cold and collisionless matter component, called dark matter (DM), which dominates over the baryon content in the matter dominated era of our Universe. The Planck satellite measurements of the Cosmic Microwave Background (CMB) temperature anisotropies have nowadays determined the amount of dark matter to an unprecedented precision, reaching the level of subpercentage accuracy in the observational determination of the abundance when combining CMB and external data Ade et al. (2016); Aghanim et al. (2018), e.g. measurements of the baryon acoustic oscillation.
Interestingly, astrophysical observation and structure formation on subgalactic scales might point towards the nature of dark matter as velocitydependent selfinteracting elementary particles. On the one hand, observations of galaxy cluster systems, where typical rotational velocities are of the order , set the most stringent bounds on the selfscattering cross section to be less than in the bullet cluster Randall et al. (2008) (in order to guarantee the production of elliptical halos MiraldaEscude (2002); Peter et al. (2013)). On the other hand, a DM selfscattering cross section of the order on dwarfgalactic scales, where velocities are of the order , would lead to a compelling solution of the cuspcore and diversity problem without strongly relying on uncertain assumptions of modelling the barionic feedbacks in simulations. This velocitydependence of the selfscattering cross section can naturally be realized in models where a light mediator acts as a longrange force between the dark matter particles. For a recent review on selfinteracting DM, see Tulin and Yu (2018).
Generically, longrange forces can lead to sizable modifications of the DM treelevel annihilation cross section in the regime where the annihilating particles are slow. For the most appealing DM candidates, known as WIMP Dark Matter Lee and Weinberg (1977); Ellis et al. (1984); Arcadi et al. (2018), such that the relic abundance in the Early Universe is set by the thermal freezeout mechanism when the DM is nonrelativistic, these effects can be sizable already at the time of chemical decoupling. Then the inclusion of the longrange force modification in the computation of the relic abundance is necessary to reach the required level of the accuracy set by the Planck precision measurement Ade et al. (2016); Aghanim et al. (2018). If the light mediators induce an attractive force between the annihilating DM particles, the total cross section is typically enhanced Hisano et al. (2004, 2007) which is often referred as the Sommerfeld enhancement Sommerfeld () or Sakharov enhancement Sakharov (1948). Additionally, such attractive forces can lead to the existence of DM boundstate solutions Pospelov and Ritz (2009); MarchRussell and West (2009); Shepherd et al. (2009). This opens the possibility for conversion processes between scattering and bound states via radiative processes, influencing the evolution of the abundance of the stable scattering states during the DM thermal history. DM scenarios with Sommerfeld enhancement or with bound states have been widely studied in the literature Belotsky et al. (2008); ArkaniHamed et al. (2009); Slatyer (2010); Feng et al. (2010); von Harling and Petraki (2014); Petraki et al. (2015); An et al. (2016); Asadi et al. (2017); Johnson et al. (2016); Hryczuk et al. (2011); Beneke et al. (2016); Freitas (2007); Hryczuk (2011); Harigaya et al. (2014); Harz et al. (2015); Beneke et al. (2015); Ellis et al. (2016); Kim and Laine (2017); El Hedri et al. (2017); Liew and Luo (2017); Mitridate et al. (2017); Harz and Petraki (2018a); Petraki et al. (2017); Harz and Petraki (2018b); Cirelli et al. (2007); MarchRussell et al. (2008); Hisano et al. (2005); Cirelli et al. (2008, 2017); Fan and Reece (2013); Bhattacherjee et al. (2014); Ibe et al. (2015); Beneke et al. (2017); Berger et al. (2008); Blum et al. (2016); Biondini and Laine (2017, 2018); Biondini (2018) and it has been shown that the main effect of such corrections is to shift in the parameter space the upper bounds on the DM mass, otherwise the theoretically predicted DM density would be too large (overclosure bound).
Classic WIMP candidates with large corrections via Sommerfeld enhancement or bound states are particles charged under the electroweak interactions, like the Wino neutralino in supersymmetric models Jungman et al. (1996) or the first KaluzaKlein excitation of the gauge boson in models with extra dimensions Servant and Tait (2003). For the supersymmetric case it was realized very early on by Hisano et al. (2004, 2007) that the Sommerfeld effect reduces the Wino density up to 30 % and pushes the mass of Wino Dark Matter candidate to few TeVs in order to obtain the correct relic density. These studies have later been extended to the case of general components of neutralino Hryczuk et al. (2011); Beneke et al. (2016). Similar and even stronger effects from the Sommerfeld enhancement and bound states were found in the case of coannihilation of the WIMP particle with charged or colored states Freitas (2007); Hryczuk (2011); Harigaya et al. (2014); Harz et al. (2015); Beneke et al. (2015); Ellis et al. (2016); Kim and Laine (2017); El Hedri et al. (2017); Liew and Luo (2017); Mitridate et al. (2017); Harz and Petraki (2018a). If the electroweakly charged Dark Matter is sufficiently heavy, the Sommerfeld enhancement or the presence of bound states due to the exchange of electroweak gauge or Higgs bosons, see e.g. Petraki et al. (2017); Harz and Petraki (2018b), are very generic as it was shown for example in the minimal Dark Matter model Cirelli et al. (2007); Mitridate et al. (2017) and in Higgs portal models MarchRussell et al. (2008). In these cases, longrange force effects play an important role also for the indirect detection limits Hisano et al. (2005); Cirelli et al. (2008); Pospelov and Ritz (2009); MarchRussell and West (2009); Cirelli et al. (2017) and especially for the Wino the Sommerfeld enhancement has lead to the exclusion of most parameter space Fan and Reece (2013); Bhattacherjee et al. (2014); Ibe et al. (2015); Beneke et al. (2017). Note that this effect can be important also when the Dark Matter is not itself a WIMP, but it is produced by WIMP decay outofequilibrium, like in the SuperWIMP mechanism Covi et al. (1999); Feng et al. (2003). Indeed in such scenario, the DM inherits part of the energy density of the mother particle and so any change in the latter freezeout density is directly transferred to the superweakly interacting DM and can relax the BBN constraints on the mother particle Berger et al. (2008).
While a lot of effort has been made to compute quantitatively the effects of a longrange force on the DM relic density employing the classical Boltzmann equation method, it is still unclear if that is a sufficient description. Indeed, considering the presence of a thermal plasma on the longrange force leads on one side to a possible screening by the presence of a thermal mass, on the other to the issue if the coherence in the (in principle infinite) ladder diagram exchanges between the two slowly moving annihilating particles is guaranteed. Moreover, from a conceptional point of view, there is yet no consistent formulation in the existing literature of how to deal with longrange forces at finite temperature, especially if the dark matter is, or, enters an outofequilibrium state (already the standard freezeout scenario is a transition from chemical equilibrium to outofchemicalequilibrium). The main concern of our work will be to clarify conceptional aspects of the derivation and the solution of the number density equation for DM particles with attractive longrange force interactions in the presence of a hot and dense plasma background, starting from first principles. From the beginning, we work in the realtime formalism, which has a smooth connection to generic outofequilibrium phenomena.
The simplified DM system we would like to describe in the presence of a thermal environment is similar to heavy quarks in a hot quark gluon plasma. For this setup it has been shown that finite temperature effects can lead to a melting of the heavyquark bound states which influences, e.g. the annihilation rate of the heavy quark pair into dilepton Burnier et al. (2008). For DM or heavy quark systems, the Sommerfeld enhancement at finite temperature has been discussed in the literature, where the chemical equilibration rate is i) estimated from linear response theory Bodeker and Laine (2013); Kim and Laine (2016, 2017) and ii) based on classical rate arguments Bodeker and Laine (2012), is then inserted into the nonlinear Boltzmann equation for the DM number density ’by hand’. Relying on those estimates, it has been shown that the overclosure bound of the DM mass can be strongly affected by the melting of bound states in a plasma environment Biondini and Laine (2017, 2018); Biondini (2018). However, strictly speaking, the linear response theory method is only applicable if the system is close to chemical equilibrium. Indeed the computation has been done in the imaginarytime formalism so far, capturing the physics of thermal equilibrium. One of our goals is to obtain a more general result directly from the realtime formalism, valid as well for systems way out of equilibrium.
Most of the necessary basics of the realtime method we use are provided in Section II as a short review of outofequilibrium Quantum Field Theory. Within this mathematical framework, an exact expression for the DM number density equation of our system is derived in Section III, where the Sommerfeldenhanced annihilation or decay rate at finite temperature can be computed from a certain component of a fourpoint correlation function. We derive the equation of motion for the fourpoint correlation function on the realtime contour in Section IV which becomes in its truncated form a BetheSalpeter type of equation. Since we close the correlation functions hierarchy by truncation, the system of coupled equations we have to solve contains only terms with DM two and fourpoint functions. In Section V, we derive a simple semianalytical solution of the fourpoint correlator under certain assumptions valid for WIMPlike freezeouts. Our result does not rely on linear response theory and it is therefore quite general applying also in case of large deviations from chemical equilibrium. The limit of vanishing finite temperature corrections is taken in Section VI, showing the consistency between our general results and the classical Boltzmann equation treatment. Here, we also compare to the linear response theory method and clarify the assumptions needed to reproduce those results. Our main numerical results for the finite temperature case are given in Section VII, both for a gauge theory and for a Yukawa potential, and discussed in detail in Section VIII. Finally, we conclude in Section IX.
Ii Realtime formalism prerequisites
The generalization of Quantum field theory on the closedtimepath (CTP) contour, or realtime formalism, is a mathematical method which allows to describe the dynamics of quantum systems outofequilibrium. Prominent applications are systems on curved spacetime and/or systems having a finite temperature. In this work, we assume that the equilibration of DM in the early Universe is a fast process, and consequently, the initial memory effects before the freezeout process can be ignored. This leads to the fact that the adiabatic assumption for such a system is an excellent approximation, motivating us to take the KeldyshSchwinger prescription^{1}^{1}1In ordinary QFT the initial vacuum state appearing in correlation functions is equivalent up to a phase to the final vacuum state. For this special situation the operators are ordered along the ’flat’ time axis ranging from to . By means of LSZ reduction formula it is then possible to relate correlation functions to the matrix and compute cross sections. This inout formalism breaks down once, for example, the initial vacuum is not equivalent to the final state vacuum. An expanding background or external sources can introduce such a time dependence. In our work, there are mainly two sources of breaking the time translation invariance. First, since we have a thermal population, we consider traces of timeordered operator products, where the trace is taken over all possible states. The many particle states are in general time dependent. Second, we have a density matrix next to the time ordering. The CTP, or, inin formalism we adopt in this work can be, pragmatically speaking, seen as just a mathematical way of how to deal with such more general expectation values. The Keldysh description of the CTP contour applies if initial correlations can be neglected and we refer for a more detailed discussion and limitation to Stefanucci and van Leeuwen (2013). of the CTP contour, as illustrated in Fig. 1. The time contour in the KeldyshSchwinger prescription consists of two branches denoted by and . The upper time contour ranges from the initial time to while the lower contour is considered to go from back to . Therefore, times on the branch are said to be always later compared to the times on . The time ordering of operator products on can generically be written as
(II.1) 
where the sum is over all the permutations of the set of operators and is the number of permutations of fermionic operators. The unit step function and the delta distribution on the KeldyshSchwinger contour is defined as
(II.2) 
Correlation functions, i.e. contour ordered operator products averaged over all states where the weight is the density matrix of the system denoted by , are defined by
(II.3) 
Let us introduce commonly used notations and properties of twopoint correlation functions of fermionic or bosonic operator pairs relevant for this work. Because of the twotime structure, there are four possibilities to align the times and on and hence four different components of a general twopoint function denoted by , where in matrix notation it can be written as
(II.4) 
Here, means and with for and the four different components of are defined as:
(II.5)  
(II.6)  
(II.7)  
(II.8) 
where on the r.h.s. of the second line applies for fermionic (bosonic) field operators. From these definitions one can recognize that not all components are independent, namely the following relation holds
(II.9) 
Furthermore, let us introduce retarded and advanced correlators defined by
(II.10)  
(II.11) 
as well as the spectral function given by:
(II.12) 
From these definitions we can derive further useful properties:
(II.13) 
In the case of free (unperturbed) propagators , the following semigroup property holds:
(II.14) 
This equality can be verified by noticing that for those time configurations the correlators are proportional to onshell plane waves in Fourier space. Note that there is no time integration in the above equation. Together with the relations in Eq. (II.13) further semigroup properties can be derived and all important ones are summarized for the use in subsequent sections in Appendix A.
As an example, the time integration over the SchwingerKeldysh contour of products of correlators can be written as:
(II.15) 
Eq. (II.15) is called a Lagereth rule and it is straightforward to work out similar rules for, e.g. different components or double integrations of higherorder products of KeldyshSchwinger correlators as they will appear later in this work. Let us move on and define Wigner coordinates according to
(II.16)  
(II.17) 
In the second line all variables are 3vectors. The Wignertransformed correlators are defined as
(II.18) 
In all computations, the tilde will be dropped such that we can write for the Fourier transformation of :
(II.19) 
One of the great advantages of separating microscopic () and macroscopic () variables according to the Wigner transformation is that Fourier transformations of integral expressions can be considerably simplified by using the gradient expansion. For example, Eq. (II.15) in Fourier space can be written as
(II.20)  
(II.21)  
(II.22) 
Throughout this work, such exponentials containing derivatives are always approximated as in the last line, defining the leading order term of the gradient expansion. For homogeneous and isotropic systems the correlators do not depend on and thus for the spatial part the leading order term is exact. For a discussion of the validity of the leading order term of the temporal part we refer to Beneke et al. (2014). Let us emphasize here, that for typical DM scenarios the leading order term is always assumed to be a very good approximation.
Next, important properties of twopoint correlators in thermal equilibrium are provided. Under this assumption, different components of correlators become related which are otherwise independent. For a system being in equilibrium (here we consider kinetic as well as chemical equilibrium), the density matrix takes the form
(II.23) 
where is the Hamiltonian of the system and factor is the inverse temperature of the system. The density matrix in thermal equilibrium can be formally seen as a time evolution operator, where the inverse temperature is regarded as an evolution in the imaginary time direction. Making use of the cyclic property of the trace it can be shown that under the assumption of equilibrium the components are related via
(II.24) 
This important property is called the KuboMartinSchwinger (KMS) relation, where () applies for twopoint correlators of fermionic (bosonic) operators. Furthermore, in equilibrium the correlators should depend only on the difference of the time variables due to time translation invariance. Consequently, the KMS condition in Fourier space reads:
(II.25) 
From this condition, various equilibrium relations follow:
(II.26)  
(II.27) 
where the FermiDirac or BoseEinstein phasespace densities are given by . Thus in equilibrium, all correlator components can be calculated from the retarded/advanced components, where the spectral function is related to those via Eq. (II.12).
General outofequilibrium observables, like the dynamic of the number density or spectral informations of the system, can be directly inferred from the equation of motions (EoM) of the corresponding correlators. Throughout this work, we derive the correlator EoM from the invariance principle of the path integral measure under infinitesimal perturbations of the fields. The equivalence of CTPcorrelators and the path integral formulation is given by
(II.28) 
and the action on the CPT contour is . collectively represents the fields, and stands for a state that could be either pure or mixed, as in Eq. (II.3). The second integral in Eq. (II.28) is a pathintegral with a boundary condition of at the initial time that we take to in the SchwingerKeldysh prescription of the contour, and the first one takes the average of with the weight of . Now, to derive the EoM for twopoint correlators, let us consider an infinitesimal perturbation , satisfying . By relying on the measureinvariance principle under this transformation, one obtains the EoM of the twopoint correlators from
(II.29)  
(II.30) 
where represents a derivative acting from the left. The same procedure can be applied to the case of , as well as for deriving the EoM of higher correlation functions. The relation between the abstract EoM of correlators and differential equations for observables will be part of the next Section. In general, a correlator EoM depends on higher and lower correlators which is called the MartinSchwinger Hierarchy. For systems where the coupling expansion is appropriate, it might be sufficient to work in the oneparticle selfenergy approximation, where the EoM are closed in terms of twopoint functions, and the kinetic equations can systematically be obtained by expanding the DM selfenergy perturbatively in the coupling constant. The kinetic equations of twopoint functions in the selfenergy approximation are also known as the KadanoffBaym equations. For example, at NLO in the selfenergy expansion of the twopoint correlators the standard Boltzmann equation is recovered.
Finite temperature corrections to nonperturbative systems, e.g. Sommerfeldenhanced DM annihilations or boundstate decays, where a subclass of higherorder diagrams becomes comparable to the LO in vacuum, are less understood in the CTPformalism. The strategy in the next section will be to address this issue by going beyond the oneparticle selfenergy approximation Beneke et al. (2014). More precisely, we derive the exact MartinSchwinger hierarchy of our particle setup in the CTPformalism by using Eq. (II.30) and truncate the hierarchy at the sixpoint function level. The system of equations is then closed with respect to two and fourpoint functions. This approach allows us to account for the resummation of the Coulomb diagrams and their finite temperature corrections at the same time. Furthermore, we show how to extract the DM number density equation from the EoM of twopoint functions and that it depends on the fourpoint correlator. The complication is that the differential equation of the fourpoint correlator is coupled to the twopoint correlator and in subsequent sections we solve this coupled systems of equations.
Iii Equations of motion in Realtime formalism for nonrelativistic particles in a relativistic plasma environment
Throughout this paper, we consider the following minimal scenario capturing important effects to study longrange force enhanced DM annihilations and bound states under the influence of a hot and dense plasma environment:
(III.1) 
The particles of the equilibrated plasma environment with temperature are the abelian mediators and the light fermionic particles with mass . Fermionic DM is assumed to be nonrelativistic, i.e. . All fermionic particles are considered to be of Dirac type. We assume the mediator to be massless, however, we provide the final results also for the case of a massive with mass .
Let us illustrate how we can get the DM nonrelativistic effective action in the thermal medium of light particles. It is obtained in two steps. First, hard modes of are integrated out. In this limit, the DM four component spinor splits into two parts, a term for the particle denoted by the twocomponent spinor and a term for the antiparticle denoted by . Second, we assume that DM does not influence the plasma environment during the freezeout process. This is typically the case since the DM number densities are smaller than that of light particles at this epoch. And thus, we may also integrate out the plasma fields by assuming they remain in thermal equilibrium. The resulting effective action on the CTPcontour for particle and antiparticle DM is given by:
(III.2) 
Dark matter longrange force interactions are encoded in the first term of the interactions in Eq. (III.2). This term includes the current and the full twopoint correlator of the electric potential on the CTPcontour which are defined by
(III.3) 
respectively. The last term in Eq. (III.2) describes the annihilation part and we only keep the wave contribution
(III.4) 
with the fine structure constant being , and summation over the spin indices are implicit. is shown in the matrix representation of the CTP formalism, see e.g. Eq. (II.4) in previous Section. Hence the delta functions on the righthandside are defined on the usual realtime axis. Similar to the vacuum theory, the annihilation part can be computed by cutting the box diagram (containing two ) and the vacuum polarization diagram (containing one and a loop of light fermions ), where now all propagators are defined on the CPTcontour. Finite temperature corrections to these hard processes in are neglected^{2}^{2}2Assuming free plasma field correlators in the computation of is a good approximation since the energy scale of the hard process is which is much larger for nonrelativistic particles than typical finite temperature corrections being of the order . Consequently, the dominant thermal corrections should be in the modification of the longrange force correlator , where the typical DM momentumexchange scale enters which is much lower compared to the annihilation scale. and for a derivation we refer to Appendix B.
In our effective action Eq. (III.2), we have discarded higher order terms in (like magnetic interactions) and also interactions with ultrasoft gauge bosons,^{3}^{3}3 To fully study the dynamics of the bound state formation Brambilla et al. (2008); von Harling and Petraki (2014); Petraki et al. (2015, 2017) at late times of the freezeout process (e.g. with being a typical binding energy), it is necessary to include emission and absorption via ultrasoft gauge bosons, e.g. via an electric dipole operator. We drop for simplicity ultrasoft contributions and discuss in detail the limitation of our approach later in this work, see Section VI.4. Note that at high enough temperature those processes are typically efficient, leading just to the ionization equilibrium among bounded and scattering states. To estimate the time when the ionization equilibrium is violated concretely, we have to take into account these processes in the thermal plasma, which will be presented elsewhere. since we focus on threshold singularities of annihilations at the leading order in the coupling and the velocity . Furthermore, our effective field theory is nonhermitian because we have integrated out (or traced out) hard and thermal degrees of freedom. The first source of nonHermite nature is the annihilation term which originates from the integration of hard degrees of freedom. A similar term would also be present in vacuum Hisano et al. (2004, 2005, 2007) and belongs to the ++ component of . Thus, as a first result we have generalized the annihilation term towards the CTP contour. Another one stems from the gauge boson propagator that encodes interactions with the thermal plasma. While the annihilation term containing in our action breaks the number conservation of DM, the interaction term containing can not. From this observation, one may anticipate that the nonhermitian potential contributions of the gauge boson propagator never lead to a violation of the DM particle or antiparticle number conservation. Later, we will show this property directly from the EoM, respecting the global symmetries of our action.
In the next sections we proceed as in the following. First, we compute the finite temperature oneloop corrections contained in the potential term explicitly. Since the number density of DM becomes Boltzmann suppressed in the nonrelativistic regime of the freezeout process, the dominant thermal loop contributions arise from the relativistic species . This implies that we can solve for independently of the DM system since we assume DM does not modify the property of the plasma. The correction terms for the DM selfinteractions are screening effects on the electric potential, as well as imaginary contributions arising from soft DM scatterings, derived and discussed in detail in Section III.1. Second, in Section III.2 the kinetic equations for the DM correlators are derived. We show how to extract from these equations the number density equation, including finite temperature corrected processes for the negative energy spectrum (boundstate decays) as well as for the positive energy contribution (Sommerfeldenhanced annihilation) in one single equation.
iii.1 Thermal corrections to potential term
In this section, we briefly summarize how the electric component of the mediator propagator gets modified by the thermal presence of ultrarelativistic fields. The plasma environment is regarded to be perturbative and in the oneparticle selfenergy framework we can write down the Dyson equation on the SchwingerKeldysh contour for the mediator:
(III.5)  
(III.6) 
where are unperturbed correlators. and are assumed to be in equilibrium and thus, according to the discussion in Section II, we only need to compute the retarded/advanced propagators. From those, we can construct all other components by using KMS condition. The Dyson equation for retarded (advanced) mediatorcorrelator can be obtained by subtracting the component of Eq. (III.5) from the component of the same equation. In momentum space this results in:
(III.7) 
where the mediator’s retarded selfenergy is defined as:
(III.8) 
A sketch of efficiently calculating the thermal oneloop Eq. (III.8) is provided in Appendix C. In the computation we utilize the Hard Thermal Loop (HTL) approximation Braaten and Pisarski (1990) to extract leading thermal corrections.^{4}^{4}4 Let us briefly summarize here the assumptions of the HTL approximation. First, we drop all vacuum contributions and only keep temperature dependent parts. Second, we assume the external energy and momentum to be much smaller than typical loop momentum which is of the order temperature (see Appendix C). The discussion of the validity of the HTL approximation depends on where the dressed mediator correlator is attached to. One can not naively argue for the case where one would attach it to the DM correlator that the external momentum is the DM momentum which is of course much larger than temperature, thus invalidating the HTL approximation by this argumentation. For example, in our case the dressed mediator correlator enters the DM singleparticle selfenergy (see Fig. 2 and Section V.3) and the dominant energy and momentum region in the loop diagram is where HTL effective theory is valid. In the HTL approximation we are allowed to resum the selfenergy contributions of the retarded/advanced component and the result is gaugeindependent. We work in the noncovariant Coulomb gauge which is known to be fine at finite temperature since Lorentz invariance is anyway broken by the plasma temperature. We find for the dressed longitudinal component of the mediator propagator in the HTL resummed selfenergy approximation and in the Coulomb gauge the following:
(III.9)  
(III.10) 
where we introduced the Debye screening mass . One can recognize that there is correction to the real part of the mediator propagator as well as a branch cut for spacelike exchange. Using the equilibrium relation Eq. (II.27), the component of in the static limit reads
(III.11) 
while for a massive mediator we have simply
(III.12) 
The static component is of special importance for describing DM longrange interactions in a plasma environment as we will see later in this work. The first term in Eq. (III.11) and Eq. (III.12) will result in a screened Yukawa potential after Fourier transformation while the second terms will lead to purely imaginary contributions. Physically, the latter part originates from the scattering of the photon with plasma fermions, leading to a damping rate Bellac (2011). Indeed in the quasiparticle picture, the mediator has a limited propagation time within the plasma, which limits as well the coherence of the mediator exchange processes. For what regards the DM particles, this term will later give rise to DM scattering with zero energy transfer, leading also to a thermal width for the DM states. In following Sections, we try to keep generality and work in most of the computations with the unspecified form and take just at the very end the static and HTL limit. Let us finally remark that the simple form of Eq. (III.12) allows us to achieve semianalytical results for the DM annihilation or decay rates in the presence of a thermal environment.
iii.2 Exact DM number density equation from correlator equation of motion
The main purpose of this section is to derive the equation for the DM number density directly from the exact EoM of our nonrelativistic action. Defining the DM particle and antiparticle correlators as
(III.13)  
(III.14) 
we derive respective EoM from the pathintegral formalism, as briefly explained at the end of Section II, for the nonrelativistic effective action given in Eq. (III.2):
(III.15)  
(III.16)  
(III.17)  
(III.18) 
The anticipated structure in Eq. (III.15)(III.18) shows the dependence of the twopoint correlators on higher correlation functions. Here, the curly brackets stand for the summation over the spin indices and is the current as already introduced in Eq. (III.3). It might be helpful to mention that we used a special property of twopoint functions of Hermitian bosonic field operators: . This exact property can be verified directly from the definition in Eq. (III.3).
In the following, the number density equation of particle and antiparticle DM is derived from this set of differential equations. First of all, we would like to clarify what is the number density in terms of fields appearing in in Eq. (III.2). For this purpose, let us switch off the annihilation term in and seek for conserved quantities. In this limit, the theory has the following global symmetries: and . The associated Noether currents for DM particle and antiparticle are:
(III.19) 
The thermalaveraged zeroth component is the number density and the average over spatial component results in the current density. We obtain the differential equation for the two DM currents directly from the twopoint function EoM, by subtracting Eqs. (III.16) from Eqs. (III.15) and Eqs. (III.18) from Eqs. (III.17), and by taking the spintrace and the limit . For the particle DM, we obtain as an intermediate result after all these steps:
(III.20) 
The trace and the curly brackets indicate the summation over the spin indices. It is important to note that the first line in the second equality cancels out, even in the case of a fully interacting correlator including finite temperature corrections. Thus, we confirm from the EoM that, e.g. nonHermitian potential corrections arising from soft thermal DM scatterings in the HTL approximation of [see Eq. (III.11)], never violate the current conservation of each individual DM species. For a homogeneous and isotropic system (vanishing divergence of current density) this would mean that the individual number densities of particles and antiparticles do not change by selfscattering processes, real physical DM scatterings, soft DM scatterings or other finite temperature corrections leading to potential contributions in .
It can be recognized, that the current conservation is only violated by the annihilation term in the last line in Eq. (III.20), since this contribution does not cancel to zero. This term can be simplified by using Eq. (III.4) and by fixing the time component to either or . We have explicitly checked that both choices of lead to the same final result. With the definition of the fourpoint correlator on the closedtimepath contour
(III.21) 
we obtain our final form of the current equations.
Current equations for particle and antiparticle :
(III.22)  
(III.23) 
The current conservation is only violated by contributions coming from . This is consistent with the expectations from the symmetry properties of the action when annihilation is turned on. Namely, only a linear combination of both global transformations leaves the action invariant which leads to the conservation of , which is nothing but the DM asymmetry current conservation. The conservation of the total DM number density, , is violated by the annihilation term.
Before we discuss the fourpoint correlator appearing in Eq. (III.22) and Eq. (III.23) in detail, let us now assume a FriedmanRobertsonWalker universe and make the connection to the Boltzmann equation for the number density that is typically adopted in the literature when calculating the thermal history of the dark matter particles. First, the spatial divergence on the left hand side of the current Eqs. (III.22) and (III.23) vanishes due to homogeneity and isotropy. Second, the adiabatic expansion of the background introduces a Hubble expansion term . Third, it can be seen from the sign of the right hand side of the current equations that only a DM loss term occurs. The production term is missing because we have assumed a priori, when deriving the nonrelativistic action, that the DM mass is much larger than the thermal plasma temperature. Within this masstoinfinity limit the DM production term is send to zero in the computation of the annihilation term and not expected to occur. Let us therefore add on the r.h.s of the current equations a posteriori the production term of the DM via the assumption of detailed balance, resulting in the more familiar number density equations:
(III.24)  
(III.25) 
The treelevel swave annihilation cross section of our system was defined as and in the last term means the evaluation at thermal equilibrium. Note that in the CTP formalism a cross section strictly speaking does not exist. The reason why this result is equal to the vacuum computation is because we computed the annihilation part at the leading order, where it is expected that zero and finite temperature results should coincide. The correlation function however is fully interacting. We summarize with two concluding remarks on our main results:

Sommerfeldenhancement factor at finite temperature: One of our findings is that the Sommerfeld factor is contained in a certain component of the interacting fourpoint correlation function, namely . This result is valid for a generic outofequilibrium state of the dark matter system. The remaining task is to find a solution for this fourpoint correlator. As we will see later, the solution can be obtained from the BetheSalpeter equation on the CTP contour, derived in the next Section. For example, expanding the BetheSalpeter equation to zeroth order in the DM selfinteractions and inserting this into Eq. (III.24) and Eq. (III.25) would result in a wellknown expression for the number density equation of the DM particles with velocityindependent annihilation. As we will see later, higher terms in the interaction or a fully nonperturbative solution contain the finite temperature corrected negative and positive energy spectrum. In other words, contains both, the bound state and scattering state contributions at the same time and at finite temperature they turn out to be not separable as it is sometimes done in vacuum computations. Bound state contributions will automatically change the cross section into a decay width and thus, appearing on the l.h.s. of Eq. (III.24) is the total number of particles including the ones in the bound states and similar interpretation for the antiparticle .

Particle numberconservation: In Section III.1, we have seen that the thermal corrections to the mediator propagator can contain, next to the real Debye mass, an imaginary contribution. It was shown that these nonhermitian corrections to the potential never violate the particle numberconservation due to the exact cancellation of the second line in Eq. (III.20). This was expected from the beginning, since, when switchingoff the annihilation , the nonrelativistic action in Eq. (III.2) has two global symmetries and . The conserved quantities are the particle and antiparticle currents in Eq. (III.22) and (III.23) in the limit (vanishing r.h.s). When annihilation is included, the nonrelativistic action is only invariant if both global transformations are performed at the same time, resulting in the conserved asymmetry current . We conclude that thermal corrections can never violate these symmetries, even not at higher loop level. On the other hand, the solution of the Sommerfeld factor is contained in and hence the annihilation rate will depend on the thermal loop corrected longerange mediator , as we will see in the next section.
Iv Twotime BetheSalpeter equations
The exact number density Eq. (III.24) depends on the Keldysh fourpoint correlation function . In this Section, we derive the system of closed equation of motion needed in order to obtain a solution for this fourpoint function, including the full resummation of Coulomb divergent ladder diagrams. The result will be a coupled set of twotime BetheSalpeter equations on the Keldysh contour as given by the end of this section, Eq. (IV.20)(IV.21). They apply in general for outofequilibrium situations and include in their nonperturbative form also the boundstate contributions if present. In order to arrive at those equations, a set of approximations and assumptions is needed. We therefore would like to start from the beginning in deriving those equations, which might lead to a better understanding of their limitation.
In the first simplification, we treat the annihilation term as a perturbation and ignore it in the following computations, since the leading order term in the annihilation part is already contained in Eq. (III.24). The exact set of EoMs for two and fourpoint functions in the limit are given by:
(IV.1)  
(IV.2)  
(IV.3) 
where we will work for the rest of this work with the conjugate antiparticle correlator and, here, the spinuncontracted correlators are defined as
(IV.4)  
(IV.5)  
(IV.6)  
(IV.7) 
The EoMs for the twopoint functions Eq. (IV.1) and (IV.2) are equivalent to Eq. (III.15) and (III.17) in the limit , and we have just rewritten them in terms of the fourpoint correlators and conjugate antiparticle propagator, defined in Eq. (IV.4)(IV.7). In our notation, the spinor indices of the operators having equal spacetime arguments in the fourpoint correlators are summed and is the current as defined in Eq. (III.3). The different conventions for the and fourpoint correlators are because and are the creation operators. From the exact differential equation of the fourpoint correlator in Eq. (IV.3), it can be seen that the correlator hierarchy is still not closed yet, since it depends on the sixpoint function. We close this hierarchy of correlators by truncating the sixpoint function at the leading order:
(IV.8) 
i.e. only the integral kernel of the six point function containing the eightpoint function was dropped. The set of correlator differential equations is now closed under this truncation procedure. A fully selfconsistent solution requires in principle to solve the equations for the five correlators