Dark Matter Sommerfeld-enhanced annihilation and Bound-state decay at finite temperature

Traditional computations of the dark matter (DM) relic abundance, for models where attractive self-interactions are mediated by light force-carriers and bound states exist, rely on the solution of a coupled system of classical on-shell Boltzmann equations. This idealized description misses important thermal effects caused by the tight coupling among force-carriers and other charged relativistic species. We develop for the first time a comprehensive ab-initio derivation for the description of DM long-range interactions in the presence of a hot and dense plasma background directly from non-equilibrium quantum field theory. Most importantly, 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 bound state decay, turns out to have in general a non-quadratic dependence on the DM number density. This generalizes the form of the conventional Lee-Weinberg equation which is typically adopted to describe the freeze-out process. We prove that our general number density equation is consistent with previous literature results under certain limits.

Traditional computations of the dark matter (DM) relic abundance, for models where attractive self-interactions are mediated by light force-carriers and bound states exist, rely on the solution of a coupled system of classical on-shell Boltzmann equations. This idealized description misses important thermal effects caused by the tight coupling among force-carriers and other charged relativistic species potentially present during the chemical decoupling process. We develop for the first time a comprehensive ab-initio derivation for the description of DM long-range interactions in the presence of a hot and dense plasma background directly from non-equilibrium 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 Sommerfeld-enhanced annihilation and bound-state 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 bound-state decay, turns out to have in general a non-quadratic dependence on the DM number density. This generalizes the form of the conventional Lee-Weinberg equation which is typically adopted to describe the freeze-out 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 on-shell Boltzmann equation treatment. So far, finite temperature corrected annihilation rates for long-range 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.

I. INTRODUCTION
The cosmological standard model successfully describes the evolution of the large-scale 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 sub-percentage accuracy in the observational determination of the abundance when combining CMB and external data [1,2], e.g. measurements of the baryon acoustic oscillation.
Interestingly, astrophysical observation and structure formation on sub-galactic scales might point towards the nature of dark matter as velocity-dependent self-interacting elementary particles. On the one hand, observations of galaxy cluster systems, where typical rotational velocities are of the order v 0 ∼ 1000 km/s, set the most stringent bounds on the self-scattering cross section to be less than σ/m DM 0.7 (0.1) cm 2 /g in the bullet cluster [3] (in order to guarantee the production of elliptical halos [4,5]). On the other hand, a DM self-scattering cross section of the order σ/m DM ∼ 1 cm 2 /g on dwarf-galactic scales, where velocities are of the order v 0 ∼ 30 − 100 km/s, would lead to a compelling solution of the cusp-core and diversity problem without strongly relying on uncertain assumptions of modelling the barionic feedbacks in simulations. This velocity-dependence of the self-scattering cross section can naturally be realized in models where a light mediator acts as a long-range force between the dark matter particles. For a recent review on self-interacting DM, see [6].
Generically, long-range forces can lead to sizable modifications of the DM tree-level annihilation cross section in the regime where the annihilating particles are slow. For the most appealing DM candidates, known as WIMP Dark Matter [7][8][9], such that the relic abundance in the Early Universe is set by the thermal freeze-out mechanism when the DM is non-relativistic, these effects can be sizable already at the time of chemical decoupling. Then the inclusion of the long-range 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 [1,2]. If the light mediators induce an attractive force between the annihilating DM particles, the total cross section is typically enhanced [10,11] which is often referred as the Sommerfeld enhancement [12] or Sakharov enhancement [13]. Additionally, such attractive forces can lead to the existence of DM bound-state solutions [14][15][16]. 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  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 [55] or the first Kaluza-Klein excitation of the gauge boson in models with extra dimensions [56]. For the supersymmetric case it was realized very early on by [10,11] 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 [26,27]. 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 [28][29][30][31][32][33][34][35][36][37][38]. 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. [39,40], are very generic as it was shown for example in the minimal Dark Matter model [37,41] and in Higgs portal models [42]. In these cases, long-range force effects play an important role also for the indirect detection limits [14,15,[43][44][45] and especially for the Wino the Sommerfeld enhancement has lead to the exclusion of most parameter space [46][47][48][49]. Note that this effect can be important also when the Dark Matter is not itself a WIMP, but it is produced by WIMP decay out-of-equilibrium, like in the SuperWIMP mechanism [57,58]. Indeed in such scenario, the DM inherits part of the energy density of the mother particle and so any change in the latter freeze-out density is directly transferred to the superweakly interacting DM and can relax the BBN constraints on the mother particle [50].
While a lot of effort has been made to compute quantitatively the effects of a long-range 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 long-range 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 long-range forces at finite temperature, especially if the dark matter is, or, enters an out-of-equilibrium state (already the standard freeze-out scenario is a transition from chemical equilibrium to out-of-chemical-equilibrium). 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 long-range force interactions in the presence of a hot and dense plasma background, starting from first principles. From the beginning, we work in the real-time formalism, which has a smooth connection to generic out-of-equilibrium 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 heavy-quark bound states which influences, e.g. the annihilation rate of the heavy quark pair into dilepton [59]. 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 [34,60,61] and ii) based on classical rate arguments [62], is then inserted into the non-linear 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 [52][53][54]. 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 imaginary-time formalism so far, capturing the physics of thermal equilibrium. One of our goals is to obtain a more general result directly from the real-time formalism, valid as well for systems way out of equilibrium.
Most of the necessary basics of the real-time method we use are provided in Section II as a short review of out-of-equilibrium 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 Sommerfeld-enhanced annihilation or decay rate at finite temperature can be computed from a certain component of a four-point correlation function. We derive the equation of motion for the four-point correlation function on the real-time contour in Section IV which becomes in its truncated form a Bethe-Salpeter 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 four-point functions. In Section V, we derive a simple semi-analytical solution of the four-point correlator under certain assumptions valid for WIMP-like freeze-outs. 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. The generalization of Quantum field theory on the closed-time-path (CTP) contour, or real-time formalism, is a mathematical method which allows to describe the dynamics of quantum systems out-of-equilibrium. Prominent applications are systems on curved space-time 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 freeze-out 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 Keldysh-Schwinger prescription1 of the CTP contour, as illustrated in Fig. 1. The time contour C in the Keldysh-Schwinger prescription consists of two branches denoted by τ + and τ − . The upper time contour τ + ranges from the initial time t i = −∞ to t = ∞ 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 C can generically be written as (1) , t p (2) ... θ C t p(n−1) , t p(n) O p(1) (t p(1) ) ... O p(n) (t p(n) ), (II. 1) where the sum is over all the permutations P of the set of operators O i and F(p) is the number of permutations of fermionic operators. The unit step function and the delta distribution on the Keldysh-Schwinger contour is defined as

II. REAL-TIME FORMALISM PREREQUISITES
Correlation functions, i.e. contour C-ordered operator products averaged over all states where the weight is the density matrix of the system denoted byρ, are defined by 1 In ordinary QFT the initial vacuum state Ω appearing in correlation functions Ω | in T [O(x, ...)] |Ω in 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 t in = −∞ to t out = ∞. By means of LSZ reduction formula it is then possible to relate correlation functions to the S matrix and compute cross sections. This in-out 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 time-ordered 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, in-in 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 [63].
Let us introduce commonly used notations and properties of two-point correlation functions of fermionic or bosonic operator pairs relevant for this work. Because of the two-time structure, there are four possibilities to align the times x 0 and y 0 on C and hence four different components of a general two-point function denoted by G(x, y), where in matrix notation it can be written as Here, G σ x σ y means x 0 ∈ τ σ x and y 0 ∈ τ σ y with σ i = ± for i = x, y and the four different components of G(x, y) are defined as: 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 Furthermore, let us introduce retarded and advanced correlators defined by as well as the spectral function given by: From these definitions we can derive further useful properties: In the case of free (unperturbed) propagators G 0 , the following semigroup property holds: (II.14) This equality can be verified by noticing that for those time configurations the correlators are proportional to on-shell 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 Schwinger-Keldysh contour C of products of correlators can be written as: 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 higher-order products of Keldysh-Schwinger correlators as they will appear later in this work. Let us move on and define Wigner coordinates according to In the second line all variables are 3-vectors. The Wigner-transformed correlators are defined as G(t, r, R, T) ≡ G(T + t/2, R + r/2, T − t/2, R − r/2) = G(t x x, t y y).
(II. 18) In all computations, the tilde will be dropped such that we can write for the Fourier transformation of G(x, y): One of the great advantages of separating microscopic (t, r) and macroscopic (T, R) 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 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 R 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 [64]. 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 two-point 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 where H is the Hamiltonian of the system and β factor is the inverse temperature T 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 This important property is called the Kubo-Martin-Schwinger (KMS) relation, where − (+) applies for two-point 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: From this condition, various equilibrium relations follow: where the Fermi-Dirac or Bose-Einstein phase-space densities are given by n F/B (ω) = 1/(e βω ± 1). Thus in equilibrium, all correlator components can be calculated from the retarded/advanced components, where the spectral function G ρ is related to those via Eq. (II.12). General out-of-equilibrium 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 CTP-correlators and the path integral formulation is given by (II. 28) and the action on the CPT contour is S = ∫ x ∈ C L(x). µ 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 path-integral with a boundary condition of µ i at the initial time t i that we take to −∞ in the Schwinger-Keldysh prescription of the contour, and the first one takes the average of µ i with the weight of ρ(µ i ). Now, to derive the EoM for two-point correlators, let us consider an infinitesimal perturbation O † = O † + , satisfying (t i ) = 0. By relying on the measure-invariance principle under this transformation, one obtains the EoM of the two-point correlators from where δ L represents a derivative acting from the left. The same procedure can be applied to the case of O, 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 Martin-Schwinger Hierarchy. For systems where the coupling expansion is appropriate, it might be sufficient to work in the one-particle self-energy approximation, where the EoM are closed in terms of two-point functions, and the kinetic equations can systematically be obtained by expanding the DM self-energy perturbatively in the coupling constant. The kinetic equations of two-point functions in the self-energy approximation are also known as the Kadanoff-Baym equations. For example, at NLO in the self-energy expansion of the two-point correlators the standard Boltzmann equation is recovered.
Finite temperature corrections to non-perturbative systems, e.g. Sommerfeld-enhanced DM annihilations or bound-state decays, where a sub-class of higher-order diagrams becomes comparable to the LO in vacuum, are less understood in the CTP-formalism. The strategy in the next section will be to address this issue by going beyond the one-particle self-energy approximation [64]. More precisely, we derive the exact Martin-Schwinger hierarchy of our particle setup in the CTP-formalism by using Eq. (II.30) and truncate the hierarchy at the six-point function level. The system of equations is then closed with respect to twoand 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 two-point functions and that it depends on the four-point correlator. The complication is that the differential equation of the four-point correlator is coupled to the two-point correlator and in subsequent sections we solve this coupled systems of equations.

III. EQUATIONS OF MOTION IN REAL-TIME FORMALISM FOR NON-RELATIVISTIC PARTICLES IN A RELATIVISTIC PLASMA ENVIRONMENT
Throughout this paper, we consider the following minimal scenario capturing important effects to study long-range force enhanced DM annihilations and bound states under the influence of a hot and dense plasma environment: The particles of the equilibrated plasma environment with temperature T are the abelian mediators A µ and the light fermionic particles ψ with mass m T. Fermionic DM χ is assumed to be nonrelativistic, i.e. M T. 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 A µ with mass m V M. 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 p M are integrated out. In this limit, the DM four component spinor χ splits into two parts, a term for the particle denoted by the two-component spinor η and a term for the anti-particle denoted by ξ. Second, we assume that DM does not influence the plasma environment during the freeze-out 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 CTP-contour for particle η and anti-particle ξ DM is given by: Dark matter long-range force interactions are encoded in the first term of the interactions in Eq. (III.2). This term includes the current and the full two-point correlator of the electric potential on the CTP-contour which are defined by respectively. The last term in Eq. (III.2) describes the annihilation part and we only keep the s-wave contribution with the fine structure constant being α i ≡ g 2 i /4π, and summation over the spin indices are implicit. Γ s 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 right-hand-side are defined on the usual real-time axis. Similar to the vacuum theory, the annihilation part Γ s can be computed by cutting the box diagram (containing two A µ ) and the vacuum polarization diagram (containing one A µ and a loop of light fermions ψ), where now all propagators are defined on the CPT-contour. Finite temperature corrections to these hard processes in Γ s are neglected2 and for a derivation we refer to Appendix B.
In our effective action Eq. (III.2), we have discarded higher order terms in ∇/M (like magnetic interactions) and also interactions with ultra-soft gauge bosons,3 since we focus on threshold singularities of annihilations at the leading order in the coupling g χ and the velocity ∇/M. Furthermore, our effective field theory is non-hermitian because we have integrated out (or traced out) hard and thermal degrees of freedom. The first source of non-Hermite nature is the annihilation term which originates from the integration of hard degrees of freedom. A similar term would also be present in vacuum [10,11,43] 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 D that encodes interactions with the thermal plasma. While the annihilation term containing Γ s in our action breaks the number conservation of DM, the interaction term containing D can not. From this observation, one may anticipate that the non-hermitian potential contributions of the gauge boson propagator never lead to a violation of the DM particle or anti-particle 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 one-loop corrections contained in the potential term D explicitly. Since the number density of DM becomes Boltzmann suppressed in the non-relativistic regime of the freeze-out process, the dominant thermal loop contributions arise from the relativistic species ψ. This implies that we can solve for D independently of the DM system since we assume DM does not modify the property of the plasma. The correction terms for the DM self-interactions 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 A. Second, in Section III B 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 (bound-state decays) as well as for the positive energy contribution (Sommerfeld-enhanced annihilation) in one single equation.

A. Thermal corrections to potential term
In this section, we briefly summarize how the electric component of the mediator propagator D gets modified by the thermal presence of ultrarelativistic ψ fields. The plasma environment is regarded to be perturbative and in the one-particle self-energy framework we can write down the Dyson equation on the Schwinger-Keldysh contour for the mediator: where S(x, y) ≡ T C ψ(x)ψ(y) 0 are unperturbed ψ correlators. A µ 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) mediator-correlator can be obtained by subtracting the +− (−+) component of Eq. (III.5) from the ++ component of the same equation. In momentum space this results in: where the mediator's retarded self-energy is defined as: Assuming free plasma field correlators in the computation of Γ s is a good approximation since the energy scale of the hard process is ∼ M which is much larger for nonrelativistic particles than typical finite temperature corrections being of the order ∼ g ψ T . Consequently, the dominant thermal corrections should be in the modification of the long-range force correlator D, where the typical DM momentum-exchange scale enters which is much lower compared to the annihilation scale. 3 To fully study the dynamics of the bound state formation [21,22,39,65] at late times of the freeze-out process (e.g. T |E B | with E B being a typical binding energy), it is necessary to include emission and absorption via ultra-soft gauge bosons, e.g. via an electric dipole operator. We drop for simplicity ultra-soft contributions and discuss in detail the limitation of our approach later in this work, see Section VI D. Note that at high enough temperature T |E B | 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.
A sketch of efficiently calculating the thermal one-loop Eq. (III.8) is provided in Appendix C. In the computation we utilize the Hard Thermal Loop (HTL) approximation [66] to extract leading thermal corrections. 4 In the HTL approximation we are allowed to resum the self-energy contributions of the retarded/advanced component and the result is gauge-independent. We work in the non-covariant 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 µ = ν = 0 of the mediator propagator in the HTL resummed self-energy approximation and in the Coulomb gauge the following: where we introduced the Debye screening mass m 2 D = g 2 ψ T 2 /3. One can recognize that there is correction to the real part of the mediator propagator as well as a branch cut for space-like exchange. Using the equilibrium relation Eq. (II.27), the ++ component of D in the static limit reads while for a massive mediator we have simply The static D ++ component is of special importance for describing DM long-range 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 [67]. Indeed in the quasi-particle 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 D(x, y) 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 semi-analytical results for the DM annihilation or decay rates in the presence of a thermal environment.

B. 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 anti-particle correlators as we derive respective EoM from the path-integral formalism, as briefly explained at the end of Section II, for the nonrelativistic effective action S NR given in Eq. (III.2): In the following, the number density equation of particle and anti-particle 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 S NR in Eq. (III.2). For this purpose, let us switch off the annihilation term Γ s → 0 in S NR and seek for conserved quantities. In this limit, the theory has the following global symmetries: η → e iθ η η and ξ → e −iθ ξ ξ. The associated Noether currents for DM particle and anti-particle are: The thermal-averaged 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 two-point function EoM, by subtracting Eqs. (III.16) from Eqs. (III.15) and Eqs. (III.18) from Eqs. (III.17), and by taking the spin-trace and the limit y → x. For the particle DM, we obtain as an intermediate result after all these steps: 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 D including finite temperature corrections. Thus, we confirm from the EoM that, e.g. non-Hermitian potential corrections arising from soft thermal DM-ψ scatterings in the HTL approximation of D [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 self-scattering processes, real physical DM-ψ scatterings, soft DM-ψ scatterings or other finite temperature corrections leading to potential contributions in D.
It can be recognized, that the current conservation is only violated by the annihilation term Γ s 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 x 0 to either τ + or τ − . We have explicitly checked that both choices of x 0 lead to the same final result. With the definition of the four-point correlator on the closed-time-path contour we obtain our final form of the current equations.
Current equations for particle η and anti-particle ξ: The current conservation is only violated by contributions coming from Γ s . 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 ∂ µ (J µ η − J µ ξ ) = 0, which is nothing but the DM asymmetry current conservation. The conservation of the total DM number density, ∂ µ (J µ η + J µ ξ ), is violated by the annihilation term.
Before we discuss the four-point correlator appearing in Eq. (III.22) and Eq. (III.23) in detail, let us now assume a Friedman-Robertson-Walker 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 H. 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 mass-to-infinity limit the DM production term is send to zero in the computation of the annihilation term Γ s 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: The tree-level s-wave annihilation cross section of our system was defined as (σv rel ) = π(α 2 χ + α χ α ψ )/M 2 and | eq 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 Γ s at the leading order, where it is expected that zero and finite temperature results should coincide. The correlation function G ++−− ηξ,s however is fully interacting. We summarize with two concluding remarks on our main results: • Sommerfeld-enhancement factor at finite temperature: One of our findings is that the Sommerfeld factor is contained in a certain component of the interacting four-point correlation function, namely G ++−− ηξ,s . This result is valid for a generic out-of-equilibrium state of the dark matter system. The remaining task is to find a solution for this four-point correlator. As we will see later, the solution can be obtained from the Bethe-Salpeter equation on the CTP contour, derived in the next Section. For example, expanding the Bethe-Salpeter equation to zeroth order in the DM self-interactions x, x) = n η n ξ and inserting this into Eq. (III.24) and Eq. (III.25) would result in a well-known expression for the number density equation of the DM particles with velocity-independent annihilation. As we will see later, higher terms in the interaction or a fully non-perturbative solution contain the finite temperature corrected negative and positive energy spectrum. In other words, G ++−− ηξ,s 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, n η 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 anti-particle ξ.
• Particle number-conservation: In Section III A, we have seen that the thermal corrections to the mediator propagator D can contain, next to the real Debye mass, an imaginary contribution. It was shown that these non-hermitian corrections to the potential never violate the particle number-conservation due to the exact cancellation of the second line in Eq. (III.20). This was expected from the beginning, since, when switching-off the annihilation Γ s → 0, the nonrelativistic action in Eq. (III.2) has two global symmetries η → e iθ η η and ξ → e −iθ ξ ξ. The conserved quantities are the particle and anti-particle currents in Eq. (III.22) and (III.23) in the limit Γ s → 0 (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 J η − J ξ . 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 G ++−− ηξ,s and hence the annihilation rate will depend on the thermal loop corrected longe-range mediator D, as we will see in the next section.

IV. TWO-TIME BETHE-SALPETER EQUATIONS
The exact number density Eq. (III.24) depends on the Keldysh four-point correlation function G ++−− ηξ,s (x, x, x, x). In this Section, we derive the system of closed equation of motion needed in order to obtain a solution for this four-point function, including the full resummation of Coulomb divergent ladder diagrams. The result will be a coupled set of two-time Bethe-Salpeter equations on the Keldysh contour as given by the end of this section, Eq. (IV.20)-(IV.21). They apply in general for out-of-equilibrium situations and include in their non-perturbative form also the bound-state 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 Γ s 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 four-point functions in the limit Γ s → 0 are given by: where we will work for the rest of this work with the conjugate anti-particle correlatorḠ ξ and, here, the spin-uncontracted correlators are defined asḠ The EoMs for the two-point functions Eq. (IV.1) and (IV.2) are equivalent to Eq. (III.15) and (III.17) in the limit Γ s → 0, and we have just rewritten them in terms of the four-point correlators and conjugate anti-particle propagator, defined in Eq. (IV.4)-(IV.7). In our notation, the spinor indices of the operators having equal space-time arguments in the four-point correlators are summed and J is the current as defined in Eq. (III.3). The different conventions for the ηη and ξξ four-point correlators are because η † and ξ are the creation operators. From the exact differential equation of the four-point correlator in Eq. (IV.3), it can be seen that the correlator hierarchy is still not closed yet, since it depends on the six-point function. We close this hierarchy of correlators by truncating the six-point function at the leading order: i.e. only the integral kernel of the six point function containing the eight-point function was dropped. The set of correlator differential equations is now closed under this truncation procedure. A fully self-consistent solution requires in principle to solve the equations for the five correlators G η , G ξ , G ηξ , G ηη , G ξ ξ simultaneously. This is beyond the scope of this work and we have to further approximate the system in order to obtain at least a simple semi-analytical solution at the end. The solution of our target component G ηξ can formally be decoupled from the solution of G ηη and G ξ ξ by approximating the latter two quantities as the leading order contribution (dropping the integral kernel, known as Hartree-Fock approximation). Then, one can recognize that In both equations, the last step is a strict equality only if the DM particle G η (z, z) = −n η and anti-particleḠ ξ (z, z) = −n ξ number densities are equal. This is true if there is no DM asymmetry which we will assume throughout this work. In the last approximation, we perform a coupling expansion of Eq. (IV.1)-(IV.3). After inserting the results of the two-point functions into the equation of G ηξ , by using the relations Eq. (IV.9) and (IV.10) in the free limit, we obtain for the four-point correlator to the leading order in g χ : where in the last two terms, particle and anti-particle are disconnected and we have introduced the single-particle self-energies according to The first integral term in Eq. (IV.11) contains the ladder diagram exchange between particle and anti-particle. Similar equations for G ηη and G ξ ξ can also be obtained by applying the same steps. In order to obtain the spectrum of bound state solutions as well as a fully non-perturbative treatment of the Sommerfeld-enhancement we have to resumm Eq. (IV.11) somehow.

We define our resummation scheme of the four-point correlator by resumming the ladder exchange, as well as the self-energy contributions in Eq. (IV.11) on an equal footing.
In other words, we are resumming the leading order terms in the coupling expansion of the four-point correlator G ηξ and similar for G ηη and G ξ ξ . On the one hand, this is a critical point, since this procedure is not exact and we cannot guarantee that other contributions do not play an important role as well, e.g. one of the limitations are systems with large coupling constants or large DM density where we can not decouple the solution of G ηξ from G ηη and G ξ ξ . The latter limitation might not be a problem since when focussing on the freeze-out the DM density becomes dilute. On the other hand, as we will discuss in detail later in this work, the final result based on this resummation scheme behaves physically, fulfils KMS condition in equilibrium5, reproduces the correct vacuum limit, enables us to study bound states and DM Sommerfeld-enhanced annihilation at finite temperature, and in some other limits we recover the literature results based on linear response theory. Furthermore, a combined resummation of one-particle self-energies and ladder-diagram exchanges seems to be necessary at finite temperature, since similar combinations would also occur when calculating the effective potential from Wilson-loop lines [68] guaranteeing the gauge invariance.
Before we can resum Eq. (IV.11), we need some further exact rearrangements and manipulations of the four-point function components. Since we are interested in the solution of G ++−− ηξ at equal times, it is sufficient to only consider two-time four-point functions, where we will adopt the short notation G ηξ (t, t ) = G ηξ (txy, t zw). Certain combinations of the components of Eq. (IV.11) turn out to be closed, e.g. let us define 6 (IV.14) In the following we will show that, by using the semigroup properties of free correlators, the following structure of the retarded equations can be achieved and similar for the advanced. The precise terms contained in the two-particle self-energy Σ ηξ will be given later. From the form of Eq. (IV.15) it is clear that our resummation scheme as described above is just the replacement of the free two-point correlators at the end by Let us in the following sketch the way how to obtain this resummable structure. The first term on the r.h.s in Eq. (IV.15) can be obtained by using simple relations 16) In the last step we used G R G A = 0 for equal time products and when multiplying this with the unit step function as in the definition Eq. (IV.13) and (IV.14) this just projects out the respective product. The integral term is more complicated. Let us consider for simplicity only the last integral term in Eq. (IV.11) and perform the sum over the different contributions of the components for the retarded two-time four-point correlator, resulting in: where the retarded one-particle self-energy is defined as Σ R = Σ ++ − Σ +− . Since all propagators are free, we can use the semigroup property (see also Appendix A) For brevity, we suppress the space integration here. This property can be used in Eq. (IV.17) since the times satisfy the inequality due to the product of retarded correlators, resulting in: Comparing with Eq. (IV.15), we indeed find the anticipated structure. Applying similar steps to all integral terms in Eq. (IV.11), we find the following Two-time Bethe-Salpeter equations7.

Two-Time Bethe-Salpeter equations:
The products of free correlators are defined as and similar for the other components, e.g. G ++−− ηξ = G +− η,0 G +− η,0 . Furthermore, we introduced the two-particle self-energies according to 7 Similar equations where also obtained in [69], although the derivation, particle content and potential is slightly different compared to our. Nevertheless, further helpful steps for bringing the integral terms into a resumable form by using semigroup properties can be found in the Appendix of [69].
where the retarded single-particle self-energy is defined in terms of the components of the definition Eq. (V.20), namely i . The advanced two-particle self-energy is given by 24) and the most important statistical components are defined by A graphical illustration of the resummation scheme for the retarded or advanced component is shown in Fig. 2. Below, we . The terms in the brackets belong to one-particle selfenergy contributions as well as the mediator exchange between particle and antiparticle. Dots represent terms containing the DM distribution in the two-particle self-energy. Due to the Boltzmann suppression at the freeze-out, those contributions will be dropped later (see DM dilute limit in Section V B).
summarize the essential properties of our two-time BS equations, based on our resummation of one-particle self-energies as well as the mediator exchanges on an equal footing.
• Two-time structure The remarkable result of our resummation scheme is that we achieved a two-time dependence of the Bethe-Salpeter Eq. (IV.20) and (IV.21) without assuming a static property of mediator correlator D or, more general, without assuming any particular form.
• KMS condition In general, the two-time correlators G ++−− ηξ (t, t ) and G −−++ ηξ (t, t ) are related in equilibrium via the KMS condition: . Now, one of the great features of our resummation scheme is that it respects this KMS condition which is non-trivial since the equations are not exact. This means that the left and right hand side of BS Eq. (IV.20) for the components G ++−− ηξ (t, t ) and G −−++ ηξ (t, t ) transform respectively into each other in equilibrium. This can be seen, by Fourier transforming the kernels, needed for proper analytic continuation of time, and by using the fact that in equilibrium all quantities only depend on the difference of the time variable. Indeed one finds that always the statistical part in the three kernels of Eq. (IV.20) transform into their counterparts. The solution of our target component G ++−− ηξ becomes a very simple expression when utilizing the power of the KMS condition as will see in the next section.
• Similar to the two-point functions, two-time retarded and advanced components of the four-point correlator are related by , as can be shown directly from the definition Eq. (IV.13) and Eq. (IV.14).
• Other components of four-point correlators not listed above, can be constructed by the given ones [69].
• The full set of equations are able to describe Bose-Einstein condensation, e.g. relevant for fermionic systems where boundstate solutions exist and the density and chemical potential are in a critical regime. Since we focus on multi-TeV particles, the density will be always low enough to ignore those quantum-statistical effects.

V. TWO-PARTICLE SPECTRUM AT FINITE TEMPERATURE
For general out of-equilibrium situations, the coupled system of two-time BS Eqs. (IV.20)-(IV.21) might require a fully numerical treatment. However, when relying on some well-motivated assumptions which are guaranteed for WIMP like freezeouts, we show in this chapter that the coupled equations can by drastically simplified. The result will be a formal solution of our target component G ++−− ηξ , appearing in our main number density Eqs. (III.24)-(III. 25), in terms of the DM two-particle spectral function. Furthermore, we fully provide the details for finding the explicit solution of the DM two-particle spectral function from a Schrödinger-like equation with an effective in-medium potential, including thermal corrections. For a better understanding of the limitation of this final solution we share, step by step, the approximations needed. All assumptions leading to this simple result will be made systematically and clearly separated section-wise in this chapter.

A. Formal solution in grand canonical ensemble
We assume the DM system to be in a grand canonical state where the density matrix takes the form: where • = η, ξ and we assume symmetric DM resulting in equal chemical potentials µ = µ η = µ ξ . Further, we assume the Hamiltonian to commute with the number operators by treating the annihilation Γ s as a perturbation. This is valid as long as the processes driving the system to a grand canonical state are much more efficient compared to annihilation or the decay of bound states. In this sense, the chemical potential is effectively time dependent. It is related to the total number density, appearing on the l.h.s in Eqs. (III.24)-(III.25). The time dependence of the number density is set by the Hubble term and the production and loss terms appearing on the r.h.s.. Let us insert the grand canonical density matrix into the components + + −− and − − ++ and show how they are related.
, one can derive the KMS relations for the two-time four-point correlators in the presence of the chemical potentials. Recalling that the Hamiltonian is a generator of the time evolution, one finds the KMS condition for a grand canonical state: Its Fourier transform reads where we introduced the Wigner-transformed four-point correlators according to Here, we have used the fact that the operatorρ commutes with the HamiltonianĤ and the translation operatorP . Defining the two-particle spectral function as our target component is formally solved in terms of the spectral function and chemical potential by utilizing the KMS relation for a grand canonical state.
Collision term for grand canonical state: In the second line we approximated the Bose-Einstein distribution f B as Maxwell-Boltzmann, assuming the DM system to be dilute T (M − µ). The same limit should also be taken in the explicit solution of the spectral function, as it is done in the next Section. In the last equality (V.8), we have used that the spectral function only depends on E ≡ ω − P 2 /4M (as we will explicitly see later) and adopted the loose notation G ρ ηξ (0, 0; ω, P ) = G ρ ηξ (0, 0; E). As a consequence of Fourier transformation, the energy integration in Eq. (V.8) ranges from minus infinity to plus infinity.
If the theory has bound-state solutions in the spectrum, the two-particle spectral function has strong contributions at particular negative values of E (binding energy). These contributions are further enhanced by the Boltzmann factor e −βE compared to the scattering solutions with positive energy, as can be directly seen from Eq. (V.8). This is expected from the assumption of a grand canonical system. All DM states of energy E must be populated with the Boltzmann factor e −βE for a given DM number density. As a result, the bound states are preferred compared to the scattering states because their energies are smaller.
Under the key assumption of the grand canonical ensemble, our main number density Eqs. (III.24)-(III.25) are formally closed. This is because the effective chemical potential appearing in Eq. (V.8) is related to the total number density by Legendre transformation. On the one hand side, the validity of adopting a grand canonical state requires scattering processes to be efficient in order to keep DM in kinetic equilibrium with the plasma particles A µ and ψ. For light mediators this is indeed guaranteed for times much later than the freeze-out. On the other hand, if the theory has bound-state solutions and is described by a grand canonical ensemble with a single chemical potential as we have introduced, there appears a hidden assumption on the internal chemical relation between scattering and bound state contributions. As we will see later, it automatically implies the Saha condition for ionization equilibrium. Ionization equilibrium can only be achieved by efficient radiative processes like the ultra-soft emissions of mediators. In the description of our theory we have traced out from the beginning these contributions but can now formally include them by assuming ionization equilibrium. Thus, a grand canonical description of systems where bound states exists is only appropriate if the ionization equilibrium can be guaranteed. In Section VI D, we come back to this issue in detail, provide an explicit expression for the chemical potential and prove the implication of ionization equilibrium.
We would like to finally remark that once a grand canonical picture is appropriate, all finite temperature corrections to the annihilation or decay rate are contained in the solution of the two-particle spectral function in Eq. (V.8) through Eqs. (III.24) and (III.25). Indeed as we will see, the negative and positive energy solutions of the spectral function will merge continuously together if finite temperature effects are strong. In this case it turns out to be impossible to distinguish bound from scattering state solutions. Due to the form of Eq. (V.8) it is, however, not required to distinguish between these two contributions. Integrating the spectral function over the whole energy range automatically takes into account all contributions. In summary, for a grand canonical ensemble, the finite temperature corrections to Sommerfeld-enhanced annihilation and bound-state decay processes are contained in the two-particle spectral function and an explicit solution of the latter quantity is derived in following sections.

B. Two-particle spectral function in DM dilute limit
A key observation is that we have factored out in (V.8) the leading order contribution in the DM phase-space density and the remaining task is to find the two-particle spectral properties. For the computation of a two-particle spectral function we can now approximate the DM system to be dilute which is the limit T (M − µ): where G ρ is the single particle spectral function as introduced in Section II. In the DM dilute limit, the two-particle spectral function is related to the the imaginary part of the dilute solution of the retarded four-point correlator, where the result is given below.
Two-particle spectral function in the DM dilute limit: Here, we defined G ret ηξ and G adv ηξ as the DM dilute limit of the Eq. (IV.21) for G R ηξ and G A ηξ , respectively. For the rest of this section, the computation of the dilute limit of these equations is given, proving the claim ]. Applying the DM dilute limit Eq. (V.9) to the two-time BS Eq. (IV.21) for G R ηξ , one finds: where Σ ret ηξ is the dilute limit of the two-particle self-energy Σ R ηξ in Eq. (IV.23). All space-time dependences remain the same as in Eq. (IV.21) but we suppress them hereafter for simplicity. Similar limit can be taken for the advanced component. Important to note is that due to the DM dilute limit, it can be recognized that the retarded Eq. (V.11) and hence the two-particle spectral function are independent of the DM number density. For the freeze-out process the DM dilute limit is an excellent approximation. Now to continue with the proof of where in the step to the last equality we used two-time BS Eq. (IV.21) backward. The statistical two-particle self-energy in the dilute limit is given by where we have used: The first similarity is a consequence of the dilute limit. In the last equality we used G R G A = 0 for equal time products. The integral term in Eq. (V.14) vanishes by noting that in the dilute limit we have indeed −Σ ret

C. Retarded equation in static potential limit
In the previous section, we related the two-particle spectrum to the solution of the retarded equation in the DM dilute limit: . As a final step, we further simplify the two-time Bethe-Salpeter Eq. (V.11) for G ret ηξ by taking the proper static limit of the mediator correlator D, resulting finally in a simple Schrödinger-like equation. We start by acting with the inverse two-particle propagator from the left on Eq. (V.11), arriving at Here, we suppress spin-indices for simplicity but quote the final result in full form later. Let us simplify the interaction kernel in Fourier space, where we take Wigner transform in time τ ≡ t − t : Taking this Fourier transform of the kernel leads to two distinct parts: where I 1 results from the sum of the two one-particle self-energy contributions, whereas I 2 originates from the exchange term between particle and antiparticle (see Fig. 2): We can perform the two ω integrations over one-particle spectral function, where in the free limit they are given by: Now the integrals are reduced to where Ω i contain the respective on-shell energies. Noticing that the Fourier transform of D ++ (t, r) = θ(t)D −+ (t, r)+θ(−t)D +− (t, r) is given by we can now take the proper static limit Ω i → 0 which results in Introducing Wigner-momenta and Fourier transforming back with respect to the difference variables we finally end-up with the Schrödinger-like equation for the retarded four-point correlator in the static limit of the potential: Now we see that the retarded and hence also the spectral function only depends on E ≡ ω − P 2 /(4M).
Retarded two-time BS equation in static and dilute limit: where the effective in-medium potential is defined as The spectral function, we would actually like to compute, is obtained from the solution of this equation according to the relation G ρ ηξ (0, 0; E) = 2 iG ret ηξ (0, 0; E) , as proven in the previous section. In the next Section, we will derive the explicit solution of the retarded equation where we will further approximate the static mediator correlator D ++ in the hard thermal loop limit, as already given in Eq. (III.11). The first term in the effective potential in Eq. (V.31) originates from the sum of the two single-particle self-energies, while the second term accounts for the mediator exchange between particle and anti-particle. The trace in the Schrödinger-like Eq. (V.30) takes the correct spin summation into account, which we have suppressed in this Section for simplicity.

D. Explicit solution in static HTL approximation
Taking the static HTL approximation of the massless mediator as derived in Eq. (III.11), the effective potential according to Eq. (V.31) results in and φ(0) = 0 and φ(∞) = 1. One can recognize the real part of the potential is corrected by the Debye mass as expected. It shifts the energy by α χ m D (twice the Salpeter correction of single particle self-energies) and screens the Coulomb potential. At the same time, the effective potential contains an imaginary part. The physical meaning of this term is scatterings of DM with light particles in the thermal plasma. If particle and anti-particle are far away, the imaginary part must be solely determined by scattering with the thermal plasma without the Yukawa force. One can also see that this is actually the case, since the imaginary part becomes twice the thermal width of single particles −iα χ T for r → ∞ (See single particle corrections in Appendix E 2, as well as Salpeter correction in Appendix E 3). This property follows from our resummation scheme, treating the DM self-energy on an equal footing with the ladder diagram exchange. Since the finite temperature corrections introduce a constant imaginary term for large distances, we can drop the i term in the following derivation of the explicit solution of Eq. (V.30). This solution will be general and contains also the correct vacuum limit, where i has to be carefully taken into account. The effective potential in Eq. (V.32) was first obtained in [68] and reproduced subsequently by other methods [65,70]. Let us remark that we derived it independently, i.e. for the first time starting from a set of two-time Bethe-Salpeter equations on the Keldysh contour.
To derive the solution, we expand G ret ηξ in terms of partial waves leading to the l = 0 (s-wave) equation: The physically relevant boundary conditions we impose on G ret ηξ,s are listed below.
• For |r − r | → ∞, G ret ηξ,s decays exponentially (as a consequence of constant imaginary potential). Note here that, since we are working in the dilute limit, the Feynman propagator and the retarded function are the same. These requirements set the form of the solution uniquely (see also [43,71]): where g >/< are the solutions of the homogeneous differential equation whose boundary conditions are given as follows: • g < (0) = 0 and g < (0) = 1.
One can explicitly check that the solution of the form (V.36) with the boundary conditions for g >/< fulfils the requirements on G ret ηξ,s as listed above.
The two-particle spectral function, needed to compute the annihilation or decay rate according to Eqs. (V.8) and (V.10), can be written in terms of the homogeneous solution as: In the last term, prime stands for the derivative with respect to r. Formally this solution is correct but might be be troublesome in the numerical evaluation, since the real part of g > (0) has a singularity due to the 1/r behavior of the effective potential. To resolve this issue, we rewrite the imaginary part of g > (0) by means of a different solution. We closely follow the discussion given in Ref. [71]. Let us define another singular solution g s whose boundary conditions are given by g s (0) = 1 and [g s (0)] = 0. The outgoing solution for r → ∞ can be expressed by a linear combination of . The fact that g > is outgoing forces it to decay for r → ∞, which determines B. And thus, one finds The physical quantity [B] we would like to compute does not depend on the choice of [g s (0)]. This is why we can get the correct result without handling the divergence of [g s (0)]. The final step is to rewrite the singular solution as .
One may also write down the expression for g > by using Eq. (V.41): . , we finally arrive at the convenient form to evaluate the imaginary part (see also [59]): (V.43) The great advantage of this general solution is that it applies for the whole two-particle energy spectrum of our theory, i.e. for negative and positive E. Here, we have introduced the tolerance δ 1 as initial value for numerical studies. Let us finally introduce dimensionless variables, expressing distances in terms of the Bohr radius x ≡ α χ Mr, and summarize the equations in terms of these units below.
S-wave part of two-particle spectral function: .
Homogeneous equations for massless and massive mediator: In practice, we use the following initial conditions which can be obtained by power series approach: where γ C applies for the Coulomb case Eq. (V.45) while γ Y can be taken for the Yukawa case Eq. (V.46). To efficiently deal with the sometimes highly oscillatory integrand in Eq. (V.44), local adaptive integrating methods are useful. In finding the initial power in x for the imaginary part we have assumed that φ(x) ∼ 1 2 x 2 for small x [59] which is only approximately true. In general, one has to carefully check if the correct value of the potential is close enough to this approximation at the initial value which we have done for the numerical results presented in subsequent Sections.

VI. DM NUMBER DENSITY EQUATION IN GRAND CANONICAL ENSEMBLE
In the previous Section, we have obtained a formal solution of the out-of-equilibrium term G ++−− ηξ,s , entering our main number density Eq. (III.24). This was achieved by assuming the DM system is in a grand canonical state, formally solving G ++−− ηξ,s in terms of the chemical potential and two-particle spectral function by KMS relation. Inserting this formal solution given in Eq. (V.8) into the main number density Eq. (III.24) results in our Master formula for the DM number density equation in a grand canonical ensemble: where a symmetric plasma 2µ η = µ η + µ ξ is assumed and (σv rel ) is the s-wave tree level annihilation cross section. The latter quantity is averaged over initial internal degrees of freedom (spin) and summed over final.
The chemical potential µ η [n η ] is a function of the total number density n η as it appears on the left hand side of our master formula. The term G ++−− ηξ,s eq is the chemical equilibrium limit µ → 0 of Eq. (V.8), given by: We presented a general method in the previous Section of how to compute the in-medium two-particle spectral function G In subsequent sections of this chapter we demonstrate how powerful our Master Eq. (VI.1) is. We self-consistently compute the component G ++−− ηξ,s (x, x, x, x) eq and the chemical potential µ. This means in both terms the same approximations should be made to obtain a well behaved number density equation.
For a better understanding, we would like to start in the next Section VI A with the simplest case of our theory by taking the zero self-interaction and zero finite temperature correction limit. This means we take g χ → 0 and g ψ → 0 in the effective in-medium potential Eq. (V.32) and compute the spectral function. Same limits are applied to the Hamiltonian entering in Eq. (VI.4) to compute the total pressure. Under these limits, our Master Eq. (VI.1) reduces to the conventional Lee-Weinberg equation, describing constant s-wave annihilation of DM.
As a next step, we allow for long-range self-interactions but neglect finite temperature corrections. This corresponds to the limit g ψ → 0, leading to the fact that the remaining term in the effective in-medium potential is the standard Coulomb or unscreened Yukawa potential. In Section VI B, the two-particle spectrum for this simple case of our theory is shown. We stress the point that only in this limit there is a direct relation between spectral function and standard expressions for the Sommerfeld-enhancement factor or the decay width of the bound states. Section VI C completes the results by computing the chemical potential for the same limit. Combining the analytic expressions for the spectral function and chemical potential, we prove our Master formula Eq. (VI.1) to be consistent with the classical on-shell Boltzmann equation treatment for vanishing thermal corrections. We also point out that adopting a grand canonical ensemble with one single time dependent chemical potential as in our master formula implies ionization equilibrium between the scattering and bound states. A detailed discussion is given on the validity of ionization equilibrium during the freeze-out process. If no bound-state solutions exist, the only limitation of our master formula is essentially kinetic equilibrium [72,73].
We relax the assumption of zero finite temperature corrections in Section VI D. This brings us to another central result of this work: a DM number density equation, generalizing the conventional Lee-Weinberg equation and classical on-shell Boltzmann equation treatment as a consequence of accounting simultaneously for DM annihilation and bound-state decay at finite temperature. However, it should be noted that this equation in Section VI D strictly speaking only applies to the narrow thermal width case and is therefore less general compared to our Master Eq. (VI.1). This means we have neglected in Section VI D imaginary-part corrections to the effective in-medium potential for the computation of the chemical potential. While we can fully account for these non-hermite corrections in the computation of G ++−− ηξ,s (x, x, x, x) eq , it remains an open question of this work of how to consistently compute the chemical potential for the broad thermal width case. The broad thermal width case for G ++−− ηξ,s (x, x, x, x) eq we compute numerically later in this work (see Section VII). Nevertheless, we demonstrate that the chemical potential and the two-particle spectral function entering the number density equation in Section VI D can be evaluated self-consistently in the narrow thermal width limit. This approach, taking leading finite temperature real-part corrections into account, is already more general of what has been computed so far in the literature. In principle, it is possible to take a nonconsistent approach and compute G ++−− ηξ,s (x, x, x, x) eq including imaginary parts in the potential while only including real-part corrections to the chemical potential. However, some care must be taken when doing so. This is because the chemical potential corrects the functional form of the number density dependence in our master equation. We discuss in more detail the possibility of taking a non-self-consistent approach by the end of Section VI D.
Finally in Section VI E, we compare our Master Eq. (VI.1) to the previous literature, relying on the method of linear response theory. Consistency is proven in the linear regime close to chemical equilibrium.

A. Recovering the Lee-Weinberg equation
We take the limit of zero self-interactions α χ → 0 while keeping the annihilation term Γ s as a perturbation. It should be emphasized again that we have to approximate the spectral function and the chemical potential both in the same limit in order to obtain a self-consistent solution. The free spectral function without self-interactions and the ideal pressure are given by: Tr [1 2×2 ]M 3/2 E 1/2 , (VI.5) Here, H 0 is the free Hamiltonian and for a derivation of this result for the two-particle spectral function directly starting from the general expression Eq. (V.44) can be found in Appendix D. The number density can be obtained from Eq. (VI.3) by using the ideal pressure: In the second equality of the first line, we find the relation between ideal number density and the non-interacting correlator G +− η,0 (x, x). The latter quantity has to be evaluated in a grand canonical ensemble, which we have done in the third equality by using KMS condition and the DM dilute limit (see Appendix E 1). The DM dilute limit should be taken in the computation of G +− η,0 (x, x) to be consistent with the computation of the spectral function. For the latter quantity we have seen in Section V B only in the DM dilute limit it is independent of the DM number density and our general solution Eq. (V.44) relies on this assumption. In the last equality of the first line we defined the conventional chemical equilibrium number density of ideal particles. Finally, we obtain from the last equality the non-interacting (ideal) chemical potential by inversion: βµ = ln[n η /n eq η,0 ]. Note that this inversion can only be done analytically if one approximates the Fermi-Dirac distribution as Maxwell-Boltzmann (which is our DM dilute limit). Entering these results of the spectral function and the chemical potential into our master formula for the DM number density Eq. (VI.1), leads to the conventional Lee-Weinberg equation for DM particles with zero self-interactions: Here, we have recovered the standard thermal averaged cross section by using the simple substitution E = Mv 2 rel /4 for the positive energy spectrum: The last equality holds for constant s-wave annihilation cross section (σv rel ) as it is the case for our model.

B. Spectral function, Sommerfeld enhancement factor and decay width for vanishing thermal corrections
We turn now to the interacting case α χ 0 and compute the two-particle spectral function. The s-wave two-particle spectral function is numerically solved according to Eq. (V.44) in the limit of vanishing finite temperature corrections and the results are shown in Fig. 3. Poles in the negative energy spectrum represent the bound states, while the spectrum is continuous for the scattering states at positive energy. In the vacuum limit one clearly sees that the scattering states can be separated from the bound state contribution at E = 0. Due to this separation, the solution of the spectral function is directly related to the Sommerfeld enhancement factor S(v rel ) and bound-state decay width Γ n . These relations are given below and now it becomes clear that the two-particle spectral function as shown in Fig. 3 is, for the vacuum case, just a convenient way of presenting all contributions simultaneously. Relation between two-particle spectral function G ρ ηξ , Sommerfeld enhancement factor S and decay width Γ n in the limit of vanishing finite temperature corrections: where E = Mv 2 rel /4 for the scattering states, E B n is the (negative) binding energy for the bound states, and (σv rel ) is the tree-level s-wave annihilation cross section. Γ n is decay width of the bound state (not to be confused with our annihilation term Γ s at the beginning of this work).
These relations can be proven directly from our general solution Eq. (V.44), see App. D for a derivation. On a first look, the spectral function in the vacuum case seems just a nice way of presentation. Instead one should emphasize that the notion of spectral function is more general and unifies the picture of scattering state annihilation and bound-state decay. This observation becomes important for the finite temperature case discussed in Section VII, where it is impossible to separate or distinguish between annihilation and decay. The spectrum includes both. Only in the absolute vacuum case a clear distinction between annihilation and decay can be made.
Coming back to Fig. 3 and now keeping in mind the relations Eqs. (VI.12)-(VI.13). There is an infinite number of exited S-bound states for the Coulomb case (left plot) with binding energy E B n = −α 2 χ M/(4n 2 ), where n is the number of the exited states and n = 1 is the ground state with lowest binding energy , shown as the pole most to the left. At small positive energies, where v rel α χ , the spectral function is constant, resulting in the familiar scaling S(v rel ) ∝ v −1 rel according to Eq. (VI.12). In the Yukawa potential case, shown in the right plot of Fig. 3, there is a finite number of bound-state solutions. For certain ratios of φ ≡ m V /(α χ M) there exist a bound-state solution with zero binding energy (E = 0). For those special cases the Sommerfeld enhancement factor scales as S(v rel ) ∝ v −2 rel for v rel m V /M α χ , called on-resonance regime, leading to an interesting observational impact on cosmology at very late times [74,75]. Roughly, those poles where the spectral function would diverge are at the on-resonance condition φ = 6/(m 2 π 2 ) where m is integer8. The on-resonance divergences give rise to partial wave unitarity violation of the total cross section, as can be seen in the right plot of Fig. 3 at E = 0. It has been pointed out in [51] that once the imaginary contribution of the annihilation part (proportional to our Γ s ) is included self-consistently in the solution of the Schrödinger-like Eq. (V.30), then the Sommerfeld enhancement starts to saturate below the unitarity limit9. This means that for some small velocity there is a transition from the divergent scaling S(v rel ) ∝ v −2 rel to S(v rel ) ∝ const. which S-wave two-particle spectral function vs. the energy E in units of typical freeze-out temperature shown for a standard Coulomb (left) and Yukawa (right) potential. The two-particle spectral function enters directly our master formula and is weighted by the Boltzmann factor for all the energy range. results in zero spectrum at zero energy. The saturation is always present if the on-resonance condition is not exactly fulfilled. The other extreme case is if φ is taken exactly in between neighbouring on-resonance values, called off-resonance. Then, S(v rel ) never scales stronger than S(v rel ) ∝ v −1 rel and at some small velocity of the order v rel m V /M the Sommerfeld enhancement factor starts to saturate and the spectral function approaches zero.

C. DM number density equation for vanishing thermal corrections
In the previous section we have proven in the limit of vanishing finite temperature corrections a relation between spectral function, standard expression of Sommerfeld enhancement factor and the bound-state decay width. Inserting these relations Eqs. (VI.12)-(VI.13) into our Master Eq. (VI.1) leads to the following differential equation for the total number density n η : (VI.14) In the limit of zero chemical potential, defining chemical equilibrium, the r.h.s. vanishes as expected. Here, we recovered the thermal averaged Somerfeld enhancement factor: The chemical equilibrium number density for the scattering states n eq η,0 was coming out as already defined in Eq. (VI.8). This outcome is fully consistent with the result one would get from integrating the Maxwell-Boltzmann equilibrium phase-space density (multiplied by spin factor 2) of non-relativistic particles. The chemical equilibrium number density of the bound-states was defined as where the mass of bound state i is M B i = 2M − |E i B | and the subscript '0' stands for ideal bound states, respectively. The term in front of the exponential in the bound state number density Eq. (VI.16) needs some further explanation. Since we have only considered s-wave contributions to the spectral function, n eq B i ,0 is the equilibrium number density of the i-th exited Para-WIMPonium. The decay, as well as the annihilation, of Ortho-WIMPonium into three A µ would be a p-wave process. To form Para-WIMPonium there is only one spin option while for Ortho-WIMPonium there are 3, consistent with the picture of having in total 4 spin degrees of freedom. Therefore, the spin factor 1 in Eq. (VI.16) comes out correctly. When carefully looking at the term in front of the exponential in Eq. (VI.16), it can be seen that the normalization of the distribution came out as like integrating the phase space density: The kinetic term P/4M of the bound state misses the correction coming from the binding energy, because the conventional normalisation would give (n eq B i ,0 ) (c) = (M B i /2π) 3/2 e −β M B i . The reason why this correction of the order O(E i B /M) does not come out as in the conventional case can be explained by how we have approximated the Martin-Schwinger hierarchy. In our paper, we expand the equation for the four-point correlator around the product of two free propagators. This is why the spectral function only depends on E = ω − P/4M as we have shown in Eq. (V.29). If we would iterate the solution, e.g. correcting the free correlators and inserting them again into the solution of the four-point correlator, we could obtain the conventional result. However, note that the correction is small for perturbative systems, since typically O(E i B /M) ∼ O(α 2 ). We are turning to the discussion of the chemical potential µ in Eq. (VI.14). The chemical potential can be obtained by inverting the number density as a function of the chemical potential. For an ideal gas description it is known that the total number density of η-particles n η , as it appears on the l.h.s. of Eq. (VI.14), is in general just given by the sum of scattering and bound state contributions: Since we have imposed a grand canonical state with only one single chemical potential, the chemical potentials for the scattering and bound states are related and therefore the number densities are not independent quantities. Assuming a grand canonical ensemble with only one time dependent chemical potential µ implies 2µ = 2µ η = µ B i , which leads to the relation The chemical potential can now be obtained from this equation by using n η,0 = n eq η,0 e βµ resulting in: Inserting this chemical potential into Eq. (VI.14) we finally end up with the Boltzmann equation for vanishing finite temperature corrections, ideal gas approximation, and the system in a grand canonical state: The equation is closed in terms of the total number density n η . Before we discuss this result in detail, let us consider the case where we would have treated all bound and scattering states to be independent. This could have been realised in Eq. (VI.14) by assigning different chemical potentials to scattering and bound states. Then we would have ended up with decoupled equations: These are the standard equations if bound and scattering states are decoupled. They might be helpful to understand Eq. (VI. 26) better. Namely, when adding the well-known Boltzmann Eqs. (VI.27) and imposing ionization equilibrium, one would end up with Eq. (VI.26). We summarize and discuss the main findings of this Section below.
• The differential Eq. (VI.26) describes the out-of-chemical equilibrium evolution of the total number density, including the reactions ηξ { AA, ψψ} (Sommerfeld-enhanced annihilation and production) and (ηξ) B { AA, ψψ} (bound-state decay and production) under the constraint of ionization equilibrium for all times. The total number density n η counts both: free particles as well as particles in the bound state. Eq. (VI.26) is equivalent to the coupled set of Boltzmann equations including soft emissions and absorptions [21] in the limit of ionization equilibrium. The equations are independent of the bound state formation or ionization cross section since the rates are, by assumption, balanced. One can also see it from a different perspective. Via Eq. (VI.26) it is very elegant to include bound state formation or dissociation processes without calculating the cross sections or solving a coupled system of differential equations. Note that when comparing our single equation to the coupled set of Boltzmann equations in Refs. [21], the last term in Eq. (VI.26) accounting for the direct production of bound states from two photon annihilation (inverse decay) was dropped. It might have little impact since bound state formation from radiative processes can be much more efficient around the freeze-out.
• Let us discuss some asymptotic regimes of Eq. (VI.26), assuming the system has bound state contributions. For this case, the ionization degree α id has a non-trivial dependence on the total number density n η , as can be seen from Eq. (VI.24). This leads to the fact that the collision term of the number density Eq. (VI.26) has in general neither a linear nor a quadratic dependence on n η , as it is for the decoupled conventional Eqs. (VI.27). However, we recover correctly the quadratic and linear form in some regimes discussed now. Close to the freeze-out, the temperature is much larger than the binding energy of bound states. As a consequence, the bound state contribution in the ionization degree can be neglected. This can be seen directly from Eq. (VI.23). In the high temperature regime x 1, leading to α id ∼ 1 (fully ionized) and the r.h.s of the number density Eq. (VI.26) is quadratic in n η . At late times, roughly when temperature becomes of the order the binding energy, the contribution of the bound states becomes strongly enhanced due to Boltzmann factor ∝ e β |E B | in K id (T). For this low temperature regime x 1, leading to α id 1/ n η K id (T) which is much smaller than unity. This behavior is expected since bound states with energy 2M − |E i B | are thermodynamically more favoured compared to the scattering states with energy 2M for T < |E i B |. Therefore, the equilibrium limit for low temperature is that most of the particles are contained in the lowest bound state (ground state). Inserting the low temperature behavior of α id 1/ n η K id (T) into Eq. (VI.26) results in the fact that at late times the r.h.s the of number density equation is linear in n η and the proportionality factor is effectively the bound-state decay rate. Since typically the decay rates are much larger compared to the Hubble rate and the number density equation is linear in n η , the total number density gets depleted exponentially fast in time. This means at late times, if ionization equilibrium is assumed, also the free DM particle number decreases simply because it must follow the exponential fast decay of the bound states in order to maintain the imposed ionization equilibrium. In summary, if bound-state solutions are present and one would integrate numerically Eq. (VI.26) until today, effectively no DM would remain as a consequence of the imposed ionization equilibrium. However, we know from the full coupled set of classical Boltzmann equations [21] that ionization equilibrium is not maintained for all times during the DM depletion phase caused by bound-state decays. This is because, once the decay rate exceeds the ionization rate, the system leaves the ionization equilibrium. The temperature dependence of the bound state formation (BSF) and ionization determines by how much the stable components are depleted during this critical epoch. Once the BSF rate drops below the cosmic expansion rate H the depletion stops.
• If we instead would have treated bound and scattering states separately, with independent chemical potentials, we would have obtained Eq. (VI.27). In this equation one can see that the bound states are independent of the scattering states. Since the differential equation of the former is linear in the number density on the r.h.s., at some point all bound states start to decay away. The coupling between bound and scattering states via radiative processes are not included in our theory and therefore do not appear. As we have learned, in our theory treating bound states as composite particles there is only one chemical potential. Therefore, by naively just giving bound and scattering states different chemical potentials would lead to the fact that we describe the system not with a grand canonical ensemble. Moreover, the KMS condition does not hold for this case. In order to be able to introduce different chemical potentials might require to rewrite our theory in terms of effective operators, creating only scattering or only bound states respectively. Then, there might be for every operator a conserved charge and one can associate individually different chemical potentials. We will see later that this more 'phenomenological' procedure is definitely not applicable for finite temperature case. There it becomes impossible to introduce such operators since the eigenvalues might be not well defined when non-hermitian thermal corrections are included.
• For going beyond the ionization equilibrium it is required to include ultra-soft terms from the beginning and re-derive the DM correlator EoM including those corrections. This we leave for future work and we restrict our equations to be valid as long as ionization equilibrium can be maintained. If no bound-state solution exist, our equations presented here only assume kinetic equilibrium [72,73]. If bound states exist, the validity of Eq. (VI.26) can be estimated. Once the decay rate exceeds the ionization rate at late times, the regime of out-of-ionization equilibrium starts (see also [21]). Until then, our description is valid. In the next Section, we will generalize this equation for the finite temperature case. We will assume that the thermal width in the effective potential is small. Already in this case, the number of bound states is dynamical since the screening as well as the constant real part term in the effective potential are temperature dependent. They lead for decreasing temperature to an abrupt occurrence of bound-state poles in the spectral function (see also Section VII). Therefore, at finite temperature the description via Eq. (VI.26) is insufficient. The more general description can be obtained from our formalism by going back to our Master equation, expressing the annihilation or decays in terms of the two-particle spectral function. The spectral function automatically decides for a given temperature which part of the spectrum contributes to the continuum and which part is bounded, as we will see later. Moreover, as we discuss in the next Section, it is required at finite temperature to include non-ideal contributions to the idealized Saha equation, leading to finite temperature corrections of the chemical potential. Thus, for a consistent treatment we need to compute both: long-range modified annihilation and the chemical potential to the same order of approximation.

D. Chemical potential in narrow thermal width approximation
In the previous Section, we established what is the outcome of our Master equation in the limit of vanishing thermal corrections. In this Section, we compute the chemical potential for finite g χ and g ψ . The only assumption we make is that the thermal width in the mediator correlator D is subleading. Only the real-part correction in Eq. (V.32) is kept. This is called the narrow width (or quasi particle) approximation. It leads to the fact that the correlator takes the simple form where V is the screened Yukawa potential. To now evaluate the pressure explicitly, we make use of the structure of the Hamiltonian for this potential: Substituting g 2 χ → λg 2 χ and taking the partial derivative of the partition function with respect to λ we arrive at the convenient form: Spin indices are summed over equal arguments and G ++−− ηη , G ++−− ξ ξ carry the same arguments as G ++−− ηξ . Applying KMS condition, and expressing G ++−− in terms of spectral function we arrive, after introducing Wigner coordinates and Fourier transformation, at: This equation tells us that the total pressure is equal to the ideal pressure, defined as contributions and the change of the scattering phase with respect to the energy for the scattering states. The equation is then known as Beth-Uhlenbeck formula. For the following discussion, it is however not required to explicitly give those expressions and therefore we just refer to the result in standard text books for non-ideal plasmas, see [76]10. Let us remark that one can also directly solve Eq. (VI.32) numerically via the solution of the two-particle spectral function. This is the power of Eq. (VI.32). One can solve self-consistently for G ++−− ηξ,s (x, x, x, x) eq and the chemical potential (see below) entering our Master equation just by evaluating the two-particle spectral function without specifying what is a bound or a scattering state. The two-particle spectral function automatically takes into account everything. This is because Eq. (VI.32) is exact in the narrow thermal width limit. We have just reduced the problem of evaluating the total pressure to the evaluation of a four-point correlation function.
For the moment, only the structure of Eq. (VI.32) is important. Namely, when differentiating the total pressure equation with respect to chemical potential we obtain, according to Eq. (VI.3), for the total DM number densities: The subscript '0' labels ideal number densities, i.e. n η,0 is a number density for a free DM without gauge interactions as in the previous Section. Superscript 'eq' stands for chemical equilibrium and a symmetric plasma is assumed µ η = µ ξ . We stress that n η,0 should not be confused with number densities of one quasi-particle excitation of DM in thermal plasma. This equation just tells us that the exact total DM number density, n η and n ξ , including corrections from thermal plasma and bound states can be simply expressed as Eqs. . In this sense, standard Boltzmann equations (see also previous Section) typically used in numerical codes solve for the DM number density non-self-consistently. This is because they miss the (small) non-ideal corrections coming from the scattering contributions in K(T) and evaluate the chemical potential in the ideal gas approximation. For a symmetric plasma n η,0 = n ξ,0 , Eq. (VI.34) is a quadratic equation in n η,0 and the solution is given by: where α gives the ratio between the ideal number density of free particles (no bound states and no interactions) and the total number density n η , including bound states, non-ideal and thermal corrections. This is the main difference compared to the ideal definition in the previous Section. The generalized Saha Eq. (VI.39) accounts for non-ideal contributions like self-interactions as well as for finite temperature corrections. Important to note is that we do not have to define what is a bound or a scattering contribution to the total number density. The spectral function in Eq. (VI.32) does the job automatically and takes all contributions into account when integrating over the whole energy range. Since n η,0 = n eq η,0 e βµ η , we can finally determine the chemical potential µ η as a function of total number density from Eq. (VI.39).
Assuming a grand canonical state for a system having bound-state solutions in the spectrum, automatically implies the Saha ionization equilibrium. Under this assumption, the chemical potential is set by βµ η = ln αn η n eq η,0 , (VI.40) and our Master formula for the total number density Eq. (VI.1) can be written in a fully closed form as: x, x, x) eq n eq η,0 (T)n eq η,0 (T) α 2 (n η K(T))n η n η − n eq η,0 (T)n eq η,0 (T) , (VI. 41) where Importantly, the obtained number density Eq. (VI.41) is in general not quadratic in n η contrary to the conventional Lee-Weinberg equation because the generalized ionization fraction, α(n η K(T)), exhibits a non-trivial dependence on n η . As we have discussed in the previous section, this equation contains all the number violating processes of DM, i.e. annihilations and bound-state decay. The process which dominates the decrease of DM number density is determined by how α evolves in time as we have discussed at the end of the previous section. The insights we gained in this Section are important for the understanding of our work. Let us put below these results more into context.
• Eq. (VI.41) describes the out-of-chemical equilibrium (finite µ) evolution of the total number density under the constraint of ionization equilibrium and is one of our main results. It is a generalization of the idealized vacuum Eq. (VI.26), accounting for i) finite temperature corrections to the annihilation/decay rates and ii) non-ideal corrections to the chemical potential. The non-ideal corrections to the chemical potential consist of finite temperature corrections as well as scattering contributions. The main advantage of the Eq. (VI.41) is that we do not have to define what is a bound or scattering state contribution since at finite temperature this is meaningless to do. All expressions needed in order to solve our generalized number density equation numerically can be obtained by evaluating the two-particle spectral function. Since one has to integrate over the whole energy spectrum, the result manifestly takes all contributions into account, without the need of differentiating between bound and scattering states. Later in Section VII, we show the results for the two-particle spectral function G ρ ηξ (0, 0; E)| l=0 in Eq. (VI.41) including finite temperature corrections. • Let us discuss a bit more in detail the generalized Saha Eq. (VI.39). For vanishing finite temperature corrections and the ideal gas limit, we have emphasized that it reduces to the standard expression Eq. (VI.24). Since (VI.39) looks the same as the one in the ideal gas limit aside from how K(T) depends on T, one may understand how α evolves from the discussion at the end of the previous Section. However, at finite temperature its is more complicated to precisely estimate since the number of bound states is dynamical. Furthermore, it is also more complicated to discuss the case when ionization equilibrium is broken at finite temperature due to this reason. Finite temperature effects might extend the period where the ionization rate is much larger compared to the decay rate, since ψ particles can efficiently destroy the bound state. Thus it might be true that the validity of our equations holds longer compared to the vacuum case. Another difficulty at finite temperature is, only in the limit of narrow thermal width it might be possible to estimate the validity of (VI.41). This is because in order to estimate the rates one has to define what is the decay width of the bound state, which becomes hard to define beyond the narrow thermal width limit. At late times of DM freeze-out, however, we naively expect that the thermal width becomes less important and one may estimate the rates by just including the real-part corrections (but still in this case the number of highly excited bound states can be dynamical).
• According to previous discussions, we would like to emphasize again that our number density Eq. (VI.41) is not applicable to the regime where bound-state decay rates exceed the ionization rates at late times causing an out-of-ionization equilibrium state (now disregarding the issue of how we can precisely estimate those rates at finite temperature). However, our more general equations Eqs. (IV.20)-(IV.21) do not assume ionization equilibrium and can be applied to any out-of-equilibrium state. The point is, since we have dropped for simplicity soft emissions from the beginning, there are no processes like BSF via the emission of a mediator relating the number of bound and free particles. Consequently, if one would solve the general Eqs. (IV.20)-(IV.21) numerically, with an initial out-of-ionization equilibrium state, the system would remain for all times in out-of-ionization equilibrium. It is required for future work to include soft emissions via e.g. an electric dipole operator in thermal plasma, to account for a correct description of the DM thermal history at late times. We will see, however, if ionization equilibrium can be guaranteed, our description accounts for sizable finite temperature corrections during the early phase of the freeze-out which can not be captured by the classical on-shell Boltzmann equation treatment as in [21]. Thus these different approaches are in some sense complementary. Furthermore, the approach in [21] uses only the main contribution of the ground state 1S, while via Eq. (VI.41) and (V.44) it is very elegant and efficient to include all (here only s-wave) bound state contributions (but under the assumption of ionization equilibrium).
• From this Section we learned that the standard Boltzmann equations at zero temperature are a non-self-consistent set of equations. They miss scattering contributions in K(T) which can now be fully accounted for. Bearing in mind that these contributions might be small, it might also be acceptable to adopt a non-self consistent solution of Eq. (VI.41). By this we mean one can in principle compute the chemical potential in the narrow width approximation, however in the computation of G ++−− ηξ,s (x, x, x, x) eq , the two-particle spectral function can be solved in its most general form including finite temperature width (as we present in Section VII). We emphasize again that a non-self-consistent solution might cause some troubles and care should be taken. In Appendix E 3, we discuss other non-self consistent computations of the chemical potential and point out their failure, especially for late times if bound-state solutions exist.

E. Comparison to linear-response-theory method
Our dark matter system is similar to heavy-quark pair annihilation in a thermal quark gluon plasma produced in heavy ion collisions. In literature, the annihilation rate of the heavy-quark pair into dileptons is estimated from linear response theory. Let us just quote their results in the following without deeply diving into the details. For comparison we translate the expressions for SU(3) to our U(1) by adjusting color and flavour factors and call the heavy (anti) quark dark matter (ξ) χ . The linearised Boltzmann equation around chemical equilibrium is given by where Γ chem is called the chemical equilibration rate. Γ chem can be extracted from linear response theory, assuming thermal equilibrium and a perturbation linear around chemical equilibrium. It is defined as [59][60][61][62]: where Ω chem. is a transport coefficient and χ η is the heavy DM number susceptibility. The transport coefficient is quoted as where we consider always only s-wave contributions and the tree level annihilation cross section is defined as in our case, while ρ ηξ is also called a spectral function but might be defined slightly differently than ours. This expression looks similar to our term in chemical equilibrium 2(σv rel )G ++−− ηξ,s (x, x, x, x)| eq = 2(σv rel ) The number susceptibility is defined as the response of the total number density with respect to infinitesimal variation of the chemical potential: (VI.47) In the dilute limit, using Eq. (VI.34), we have χ η = βn eq η . So we find in total for the linearised equation using linear response theory: In the following, we compare this equation based on linear response theory with our Eq. (VI.41). We start from our more general expression and reproduce Eq. (VI.48) in the linear regime around chemical equilibrium while clarifying the underlying approximations. For this purpose, we have to linearise our equation around n η ∼ n eq η , where 'eq' labels chemical equilibrium. On the one hand, as we have shown in Sec. V, the spectral function, G ρ ηξ , does not depend on n η as long as the DM number densities are dilute. On the other hand, the generalized ionization fraction, α, does depend on n η and hence it must also be approximated as to be close to chemical equilibrium. One may easily see this by looking at the definition of α given in Eq. (VI.39). From this expression it is clear that if the DM number densities are close to chemical equilibrium, the ionization fraction is always close to one, i.e. α eq 1, because n eq η,0 ∝ e −β M and K(T) ∝ e β |E B | . Hence, in the linear regime, the total DM number densities are fully ionized and can be approximated as the free scattering contributions: n eq η n eq η,0 . From these observations, we reproduce Eq. (VI.48) from our Eq. (VI.41) in the linear regime around chemical equilibrium.
The linear response theory method applies, by definition, only to the linear regime around chemical equilibrium where n η ∼ n eq η . It is not possible via this method to see what is the correct form of the underlying Lee-Weinberg equation describing the non-linear regime where n η n eq η . Nevertheless, one might be tempted to use it by replacing the right-hand-side of Eq. (VI.43) with Γ chem 2n eq η (n 2 η − (n eq η ) 2 ) as done in [34,[52][53][54]61]. Most important to note is that the collision term obtained from this replacement matches our equation only if the generalized ionization fraction is close to chemical equilibrium, i.e. α eq 1, which is no longer true at late times. We can see how α depends on time from Eq. (VI.39). As we have already discussed in detail at the end of Secs. VI C and VI D, at late times of DM freeze-out the total comoving number density approaches a constant value while K(T) starts to grow rapidly for T < |E B |. As a result, we find the generalized ionization fraction to be α 1/ n η K(T) 1 at late times, which invalidates the naive replacement above. One can understand this behavior intuitively because the bound states of energy 2M − |E B | are thermodynamically favored compared to the scattering states of energy 2M at T < |E B |. For a fixed total number of DM, the bound states dominate over the scattering states at some point as in the case of the recombination. There is also another issue when correcting n eq η in Γ chem 2n eq η (n 2 η − (n eq η ) 2 ) only by the Salpeter term. Since this discussion requires some detailed knowledge about thermal corrections, we share it in Appendix E 3. In summary, we would like to emphasize that care must be taken if one matches the equation obtained by linear response theory to non-linear differential equations.
Instead, one may match the equation obtained by linear response theory to our corrected form of the Lee-Weinberg equation in the non-linear regime. By this procedure it is now also clear what the limitation exactly is. As we have discussed in detail in the previous section, our master formula for the number density equation is valid as long as ionization equilibrium can be maintained. Ionization equilibrium is broken if the decay rate exceeds the ionization rate. The temperature where this happens can be estimated for finite temperature systems, at least in the narrow width case.

VII. NUMERICAL RESULTS FOR TWO-PARTICLE SPECTRAL FUNCTION AT FINITE TEMPERATURE
We turn to the numerical solution of the two-particle spectral function for the full in-medium potential. The effects of the finite temperature corrections can be simplest understood for the case of the Coulomb potential as given in Eq. (V.32). The first correction is a real constant term that shifts effectively only the energy by α χ m D . When only taking this correction into account one would thus expect that the infinite number of bound states in the spectral function of the Coulomb case just move to lower binding energies and similar shift to the threshold as well as to the positive energy spectrum. The second real-part correction is an exponential screening of the Coulomb potential with radius m D . This introduces another effect. It leads to a disappearance of the bound states closest to the threshold since Yukawa potentials have only a finite number of bound states. The disappearance of bound states wins against the move of the poles towards lower energies at the Mott transition, where all bound states disappear and the spectrum is exclusively continuous. Additionally, we have imaginary-part corrections coming from the soft DM-ψ scatterings, leading to a finite thermal width of the bound states. Once the thermal width is comparable to the binding energy, the bound-state poles are strongly broadened.
The combination of all effects are shown in Fig. 4, where we present the numerical solution of the two-particle spectral function for the full in-medium potential according to Eq. (V.44). We show the case of the Coulomb and Yukawa potential. In this figure, we have fixed the temperature to T = M/30 (slightly below the typical DM freeze-out temperature) and varied the Debye mass where the maximal value shown corresponds to the equal charge case g ψ = g χ of our minimal model: m 2 D = g 2 χ T 2 /3. The mass of the DM is fixed to 5 TeV and the coupling α χ = 0.1 very roughly chosen to account for the correct order of the abundance.
All finite temperature effects together lead to a continuous melting of the bound-state poles. As can be seen, the melting of the bound states leads to the fact that even at negative energies, the spectrum is continuous at finite temperature. The reshuffling of the spectrum towards lower energies affects the rates exponentially according to Eq. (V.8). This is because the integrand (the spectral function) has more support at negative energies which is, due to Boltzmann factor, exponentially preferred in kinetic equilibrium. The whole integral stays convergent since at very low energies the spectral function starts to decrease faster than the Boltzmann factor increases. It now becomes clear that the notion of spectral function is more general compared to the vacuum case where one could separate the spectrum for bound and scattering states. Here, it is evident that such a distinction is impossible. It is also not necessary to do so since the integration of spectral function times Boltzmann factor takes already all contributions into account.
Nevertheless, let us discuss how only the positive energy spectrum is affected. According to the theorem of Levinson, the scattering phases (and hence the wave function at the origin) depend on the amount or properties of the bound states. This means that thermal modifications of the bound states automatically affect also the positive energy spectrum. The impact on the positive energy spectrum depends on the melting status of the bound states. In general, there can be both, a suppression or further enhancement of the positive energy spectrum as can be seen by carefully looking at the value around E = 0 in Fig. 4. We would like to stress that a suppression of the positive spectrum does not imply that the total rate is less. For the computation of the rate Dashed and solid lines correspond to the vacuum limit and the full in-medium potential, respectively. Effects at different temperatures are compared. Left side for typical freeze-out temperature and right plot at a typical temperature where the annihilation rate would be much smaller compared to the Hubble rate. one has to integrate the spectral function over the whole energy range, where G ++−− ηξ,s (x, x, x, x) eq becomes due to the reshuffling towards lower energies exponentially enhanced.
While for the Coulomb case, the impact of the melting on the positive energy spectrum is only very little (which does not mean that the overall effect is small), the impact for the Yukawa potential case can be much larger. In Fig. 5, we compare the positive energy solution of the Yukawa spectrum at zero and finite temperature, as a function of the mediator mass m V . The vacuum line (dashed) is obtained by solving the spectral function in the limit of vanishing finite temperature corrections. The Sommerfeld factor for this case can be obtained from Eq. (VI.12) according to: In the second line we used E = Mv 2 rel /4 for on-shell particles. At finite temperature this kinetic energy relation does not hold. We therefore use the second Eq. (VII.2) to define the Sommerfeld enhancement factor at finite temperature as shown in Fig. 5. The enhancement or suppression of the Sommerfeld enhancement factor due to the thermal effects is largest if the ground state is close to the threshold of E = 0 (around the first peak from the right, here φ ∼ 0.6). For our minimal model, it is also shown in the right plot of Fig. 5 that the whole temperature range of a typical freeze-out process can be affected. In the limit m V → 0 the Coulomb limit is recovered. Again, this does not mean that it is sufficient to just take the standard expression of the Sommerfeld enhancement factor of the Coulomb potential to describe the DM freeze-out. There is also a contribution from the negative energy spectrum. Therefore, one has to be careful in interpreting solution (Sommerfeld enhancement factor) can be suppressed or equal compared to the vacuum case, but on the other side the total G ++−− ηξ,s (x, x, x, x) eq entering our master formula Eq. (VI.1), which requires the integration over the whole energy spectrum, can be enhanced.
As an extreme example of this situation, let us discuss the case where there are no bound states (e.g. Yukawa potential in the Born regime φ 1). The finite temperature spectral function for this case is shown in the right plot of Fig. 6. Indeed, the positive energy spectrum is dominantly suppressed compared to the vacuum case but when integrating the spectral function over the Boltzmann factor in the whole energy range there is still an enhancement of order 1% compared to the vacuum case. Another extreme example, where the corrections to the positive spectrum are strongest, is the case where only the ground state exists and is close to threshold [34,52,53]. This example is shown in the left plot of Fig. 6. We find in this case, the value of the integration of the spectral function times the Boltzmann factor is by up to 10% (30%) larger compared to the vacuum case without bound-state peak at the typical freeze-out temperature T = M/30 (T = M/90). The correction increases for lower temperature due to the Boltzmann factor.
The total rate is only proportional to the integration of the spectral function times the Boltzmann factor. There are additional finite temperature corrections to the chemical potential (see previous Section) which can be also obtained from the spectral function. We have not explicitly computed those non-ideal corrections yet but leave it for future work once we have included ultra-soft emissions in our system description.

VIII. DISCUSSION
A self-interacting DM system, where long-range forces and bound-state solutions exist, is in general a complex ensemble where many processes with different rates are taking place at the same time during the DM thermal history. Essentially, there are three quite different approaches in the literature with distinct motivation to describe the evolution of the abundance of the stable components for such systems: 1. The first approach is based on a coupled set of classical on-shell Boltzmann equations. If bound-state solutions are absent, the description of the DM freeze-out acquires dominantly corrections from the Sommerfeld-enhanced annihilation of free DM particles [10,11]. If the two-particle spectrum has support at negative energies, the free DM particles can form a bound state via radiative processes [21,77]. The reverse process can also happen, called ionization. If there are several bound-state solutions present, further processes like excitation or de-excitation can happen [22,39]. All those processes are in general coupled, and as we see, the list of Boltzmann equations needed to describe such systems can be long. When relying on those classical Boltzmann equation computations, treating e.g. the number density of free particles and bound states separately and as idealized, potential strong modifications arising from higher-order plasma interactions might be missed. In this approach, however, it is always guaranteed that the non-linearity of out-of-chemical equilibrium reactions are accurately described. And there can be in general many such out-of-equilibrium reactions as listed above.
2. The second approach starts from the EoM of correlation functions on the Keldysh contour and takes into account some finite temperature corrections. The major difference to our work is that in [64] it is assumed that the correlator hierarchy can by truncated at the lowest order, resulting in closed equations for the two-point functions in terms of the one-particle self-energy only. One of the equations are the so called kinetic equations, being the differential equations for describing the evolution of observables in terms of the macroscopic Wigner coordinates. In the one-particle self-energy approximation they are also known as Kadanoff-Baym equations. Expanding the self-energy in terms of the coupling to NLO results in the standard Boltzmann equation. At NNLO first finite temperature corrections enter. The advantage of a fixed order calculation is that infrared divergences, arising at NNLO cancel [64]. At NNLO in the self-energy expansion of the kinetic equations, the thermal corrections turn out to be strongly suppressed, i.e. to high power in T/M, compared to the NLO result. One should, however, keep in mind that there are next to the kinetic equations also the equations for the microscopic Wigner coordinates, called mass-shell equations accounting for, e.g., thermal corrections to the dispersion relation. Kinetic and mass-shell equations are in general coupled. Therefore, a self-consistent solution in principle requires to take account of corrections also from the mass-shell equation. In any case, the problem within this systematic approach is that a fixed order calculation can never account for correctly describing the Sommerfeld enhancement beyond the Born regime and also bound-state solutions will never appear.
3. The third approach addresses the description of long-range force systems at finite temperature in a non-perturbative sense, i.e. by resummation of the Coulomb divergent ladder diagrams including thermal corrections. Clearly, first attempts were made in the literature of heavy quark pair annihilation in a quark gluon plasma [59], produced in heavy-ion collisions at the LHC. More recently, some of these authors have applied the same techniques also to the DM freeze-out [34]. The method is based on linear response theory [60][61][62], estimating the DM Sommerfeld-enhanced annihilation and bound-state decay from a spectral function including finite temperature corrections. It has been shown in Ref. [52][53][54] that the DM overclosure bound, computed by this method, can be strongly affected by finite temperature effects if bound-state solutions exist. Compared to a fixed order calculation as in approach 2, the finite temperature corrections are larger. The reason is because the mass-shell equations are solved by resummation of the Hard thermal loop contribution. Albeit there are potentially strong effects, the linear response theory is strictly speaking valid only for systems close to thermal equilibrium, e.g. n ∼ n eq . At finite temperature the spectral function can in general depend on the DM density. Therefore, it is a priory not clear if the transport coefficients extracted from linear response theory can be inserted into a non-linear Boltzmann equation describing the DM freeze-out in a non-linear regime where n n eq . From vacuum computations it is known that the Sommerfeld effect can still be efficient in such a regime. This is because the transition n ∼ n eq to n n eq happens in a short time, since n eq ∝ e −M/T decreases rapidly. To the best of our knowledge this method, inserting transport coefficients obtained from linear response theory into a non-linear Boltzmann equation, has not been tested so far by using other treatments applying for generic out-of equilibrium situations.
Our formalism, presented in this work, aims towards a first step in unifying the approach 1 and 3, by generalizing the approach 2 for long-range force systems. In other words, we derived from the EoM of Keldysh correlation functions the number density equation for DM including finite temperature corrections and accounting for the full resummation of Coulomb divergent ladder diagrams. This allows to study the finite temperature corrected Sommerfeld-enhanced annihilation as well as bound-state decay. Moreover, our Master Eq. (VI.1) is able to describe the correct non-linear transition to out-of-chemical equilibrium, i.e. the freeze-out process.
Although we have derived all equations on the Keldysh contour, and therefore they should be valid for any out-of-equilibrium situation of the system, the reader should be reminded what precisely our system is. While it remains true that we can describe correctly Sommerfeld-enhanced annihilation and bound-state decay in the presence of a relativistic plasma background for out-of equilibrium situations, we have dropped from the beginning, when deriving our nonrelativistic effective action, ultra-soft contributions of the fully relativistic action. Hence, Eqs. (IV.20)-(IV.21) are missing ultra-soft contributions leading to bound state formation and ionization processes via the emission or absorption of a mediator, as well as contributions to excitation or de-excitation processes if multiple bound states exist. Once ultra-soft terms are included in the system of equations, we expect the final equations, if finite temperature effects are neglected, to coincide with the full set of equations of approach 1. Moreover, the inclusion of emission and absorption in the Keldysh formalism might lead to new insights in the production rate of dileptons or photons, produced from heavy-quark pair-annihilation in a quark gluon plasma.
In the second half of the work, we have indirectly included all bound-state formation, ionization, excitation and de-excitation processes. This was achieved by assuming our system is in a grand canonical state with one single time dependent chemical potential as in our Master Eq. (VI.1). The important observation that adopting a grand canonical picture automatically implies ionization equilibrium if bound states are present was by far not obvious to us. This key observation brought us to the conclusion that our equations in the limit of vanishing thermal corrections are equivalent to the coupled system of classical Boltzmann equations in the limit of ionization equilibrium. Thus we have shown that under certain assumptions our approach and approach 1 consistently fall together. Another important point based on this observation was that, since the ionization fraction at chemical equilibrium is close to unity, our and approach 3 are equivalent to approach 1 in the regime linear near chemical equilibrium.
Important to recognize was that our approach and approach 3 give different results if the transport coefficients extracted from linear response theory is inserted into a non-linear Boltzmann equation just by replacing Γ chem (n η − n eq η ) with Γ chem 2n eq η (n 2 η − (n eq η ) 2 ) [52][53][54]. This is because the ionization fraction depends on n η where another non-linearity comes in, and in particular, the ionization fraction will be much smaller than unity at late times. This is intuitively because the bound states are exponentially favored compared to the scattering states for T < |E B |. Furthermore, the ionization fraction counteracts against the exponential grow of Γ chem or of our G ++−− ηξ,s (x, x, x, x) eq for late limes if bound-state solutions exist. In a word, while the spectral function is identical between the linear response and ours in the DM dilute limit, the ionization fraction makes the difference. This effect is non-negligible when at late times the DM gets depleted by bound state formation effects.
Our Master Eq. (VI.1) can not be used at very late times where ionization equilibrium is not maintained. Therefore, one has to be careful in relying on our so far simplified treatment for all times during the DM thermal history. More generally, when using our equations, it has to be ensured that the rates driving the system to kinetic and ionization equilibrium are much faster than any other rates leading to a potential out-of-kinetic or -ionization equilibrium state. In the case of no ψ-particles (no finite temperature corrections), it was shown in Refs. [21,22,39] that the decay of the bound state becomes faster than the ionization via emission and absorption processes by an electric dipole operator at some point, which breaks the ionization equilibrium at late times. Later, when the bound state formation becomes inefficient compared to the cosmic expansion, the dark matter number freezes out completely. Estimating the valid regime of our approach in the presence of finite temperature corrections is a more complicated task. To draw a definitive conclusion in our case, one has to estimate these processes, including emission and absorption of ultra-soft gauge bosons, in the presence of the thermal plasma. As we have discussed already in detail, this might only be realizable if the thermal width is negligible compared to the real-part corrections. Furthermore, one has to keep in mind that the number of existing bound-state solutions is temperature dependent when already only real-part corrections are taken into account.
After this warning, we now would like to discuss the case where a grand canonical description with one single chemical potential is justified. As we have in detail presented in this work, all finite temperature corrections are then totally encoded in the solution of the DM two-particle spectral function. In the presence of ultra-relativistic fermionic particles in the background, the hard thermal loop resummed corrections to the DM system can be classified into three contributions [see effective in-medium potential Eq. (V. 32)]. The first two contributions are real-part corrections to the DM effective in-medium potential. The first one leads to an energy shift (Salpeter correction) in the DM two-particle spectral function by an amount of α χ m D towards lower energies compared to the vacuum case. Second, the Debye mass m D leads to a screening of the Coulomb potential, resulting in a temperature dependent Yukawa potential with screening radius m D . The third contribution to the effective potential is purely imaginary and originates from soft DM scattering processes with the ultra-relativistic fermions in the hot and dense background.
Let us first discuss the two real-part corrections and their implication on the two-particle spectral function for the case of a Coulomb potential. In the limit of vanishing real-part corrections, it is well known that a Coulomb system has infinitely many bound-state solutions. Furthermore, the bound-state solutions and the scattering states can clearly be separated sharply at the energy E = 0 in the two-particle spectral function (see Fig. 3). If the real-part finite-temperature corrections are included, bound states close to the threshold E = 0 disappear in the spectrum. This is simple to understand. First of all, by the real-part corrections the Coulomb potential transforms into a Yukawa potential. Yukawa potentials have only a finite number of bound states. Secondly, due to the energy shift caused by the other real-part correction, the threshold is lowered. The combination of these two effects causes the disappearance of highly excited bound states close to the threshold and the spectrum is continuous instead of discrete. The effect gets stronger with increasing Debye mass leading at the Mott transition to a total disappearance of all bound states.
Already when only real-part corrections to the in-medium potential are included, the number of bound states as well as their binding energies are temperature dependent according to the discussion above. It implies that a sharp definition of bound and scattering states can not be made for all times during the DM thermal history. However, for our total number density equation, it is NOT required at all to distinguish between bound and scattering states. For example, the computation of the ionization fraction via the generalized Saha Eq. (VI.39) can always be performed without specifying what is a bound or scattering contribution. Another example is the spectral function entering the total number density in the production term. Also here, the integration of the spectral function automatically takes into account all contributions from the spectrum. Only in the absolute vacuum limit, it is possible to separate contributions. At finite temperature everything is mixed into one single object. A separation would only cause problems like unphysical jumps in projected thermodynamical quantities when bound states abruptly disappear (if only real-part corrections are included).
The mixing between scattering and bound-state solutions becomes even stronger once the imaginary contributions to the effective in-medium potential are included. As we have seen in Section VII, these corrections lead to a thermal width of the peaks of the bound states and a continuous melting of the poles for increasing Debye mass, as illustrated in Fig. 4. Note that instead of changing the independent coupling in the Debye mass, we could have also increased the number of generations of ultra-relativistic particles in the plasma. The Debye mass is proportional to the square root of number of generations, and thus, if most of the particle content of the Standard Model would run in the thermal loop, we can have a large Debye mass although the coupling is still small. The broadening of the peak and the shift towards lower binding energies increases the annihilation or decay rate exponentially (again its not necessary to distinguish between these two) due to the integration of the product of two-particle spectral function and the Boltzmann factor, as in Eq. (V.8).
Although, we focused on a simple U(1) like theory and s-wave contributions in the present work, most of the equations for higher gauge theories, scalar mediators, or higher partial waves will change in the expected way. Let us already mention some major changes. The mediator self-energy would acquire further contributions from self-interactions as well as different colour or flavour pre-factors. The definition of the effective in-medium potential in Eq. (V.31) remains the same and is computed from the specific dressed mediator correlator. The r.h.s of the Schrödinger-like Eq. (V.30) will be proportional to the number of colours. Our number density Eq. (III.24) in the case of velocity dependent tree-level annihilation cross sections (like in p-wave case) will have a space derivative on the r.h.s.. As expected, our formalism breaks down for temperatures around the confining scale of confining theories.

IX. SUMMARY AND CONCLUSION
Traditional computations of DM Sommerfeld-enhanced annihilation and bound-state decay rates rely on the assumption that reactions of such processes are taking place under perfect vacuum conditions. In this work we developed a comprehensive derivation of a more general description, taking into account non-ideal contributions arising from simultaneous interactions with the hot and dense plasma environment in the early Universe. We have derived the evolution equation for the DM number density which is applicable to the case where scattering and bound states get strongly mixed due to the influence of the thermal plasma surrounding. Our master Eq. (VI.1) for the total DM density simultaneously accounts for annihilation and bound-state decay and hence its collision term is in general not quadratic in the DM number density. We showed that finite temperature effects can lead to strong modifications of the shape of the two-particle spectrum, which in turn modifies the DM annihilation or decay rates.
The Keldysh formalism we adopted throughout this work applies for the description of the dynamics of generic out-ofequilibrium states. Within this mathematical framework, we derived in the first part of this work directly from our nonrelativistic effective action the exact equation of motion of the DM two-point correlation functions. We extracted for the first time from those EoM the differential equation for the DM number density [see Eq. (III.24) and (III. 25)], which turns out to only depend on a special component of the DM four-point function on the Keldysh contour, namely G ++−− ηξ . Let us emphasize again that this equation for the number density is exact within our nonrelativistic effective action, however not closed since it depends on the solution of this four-point correlation function. The long-range force enhanced annihilations, the decay of bounded particles as well as the finite temperature corrections are all contained in the solution of this one single four-point correlator.
In the second part of this work, we derived the EoM for the DM four-point function on the Keldysh contour. We developed the approximations needed in order to close the hierarchy of correlators but at the same time keep the resummation of Coulomb divergent ladder diagrams as well as the finite temperature corrections. Based on our approximation and resummation scheme, the final form of the equation for our target component G ++−− ηξ is physically sound and maintains important relations like the KMS condition in equilibrium. The coupled system of equations is general enough to apply for the description of DM out-of-chemical equilibrium states.
In the third part, we explored further approximations needed in order to obtain a simple solution to our target component and to reproduce from our general equations the results in the literature, based on different assumptions. So far existing literature has estimated transport coefficients from linear response theory and entered those into a non-linear Boltzmann equation by classical rate arguments [34,[52][53][54][59][60][61][62]. We have proven that our master Eq. (VI.1) is equivalent to the method of linear response only in the linear regime close to chemical equilibrium. Finally, we must point out that the Lee-Weinberg equation, adopted in [34,[52][53][54]61] to re-derive the DM overclosure bound in the non-linear regime, is not the correct form of the number density equation to use if bound-state solutions exist in the spectrum. The ionisation fraction causes the difference as discussed in great detail in our work.
When taking the vacuum limit, our master equation reduces correctly to the coupled system of classical Boltzmann equations for ideal number densities of bound and scattering states in the limit of ionization equilibrium. In our method, it came out as a consequence of assuming the system is in a grand canonical state. Namely, we have proven that the assumption of a grand canonical state automatically implies the Saha ionization equilibrium if bound-state solution exist. One has to take the assumption of ionization equilibrium to be fulfilled for all times with a grain of salt for the following reason. From the vacuum treatment it is known that the duration of Saha ionization equilibrium is limited. Therefore, when using our Master equation one has to carefully check that this condition is satisfied for a sufficiently long period. And especially when the assumption of ionization equilibrium is not justified, one has to make sure that at least the abundance of the stable scattering states are not affected by out-of-ionization equilibrium effects which might be model dependent.
The reason why in our Keldysh formalism we can not resolve this issue at the moment lies in one particular approximation, made from the beginning. Ultra-soft emissions and absorptions were dropped for simplicity. We leave the inclusion of those quantities for future work, but expect once they are included we can fully recover the general set of coupled classical Boltzmann equations in the vacuum limit of our (future) updated equations. Moreover, this would allow us to describe Sommerfeld-enhanced annihilation and bound-state decay at finite temperature for the first time beyond the ionization equilibrium.
In the regime where ionization equilibrium is maintained, we have shown that finite temperature effects strongly mix bound and scattering states and the effects are all encoded in the solution of the two-particle spectral function. Let us remark that the numerical results for the spectral function obtained in Section VII are compatible with the linear response theory approach [34,[52][53][54][59][60][61][62], although we started from a completely different method. The component G ++−− ηξ,s eq in our master Eq. (VI.1) can be enhanced by much more than 10%. In addition, our master equation is applicable to the non-linear regime beyond the limitation of linear response if, at least, the ionization equilibrium is maintained. These results make it definitely worthwhile to further generalize our Keldysh description in order to correctly describe the out-of-ionization equilibrium transition at late times by including contributions from the ultra-soft scale.

Appendix B: Annihilation term
Although one can directly compute Γ s (x, y) defined on the closed-time-path contour, it is instructive to see how one can recover it from the more common computation of annihilations. Usually, we compute the matrix element by means of Feynman correlators so as to evaluate annihilations, which means that both x and y are on the C + contour. And thus, it yields the upper left component of Γ s , i.e. Γ ++ s . The question is how to recover the remaining three components. To answer it, let us go back one step further, namely before integrating out hard products of the annihilation. Suppose that the interaction with them takes the following form; By definition, incoming energy/momentum from O s is much smaller than M, which justifies the following approximation: The inverse relations are given by: By using these relations, the one-loop expression of the retarded mediator correlator can be simplified as Tr[γ µ S ++ (x, y)γ ν S ++ (y, x)] − Tr[γ µ S +− (x, y)γ ν S −+ (y, x)] = 1 2 Tr[γ µ S R (x, y)γ ν S S (y, x)] + 1 2 Tr[γ µ S S (x, y)γ ν S A (y, x)]. (C.5) Let us take the limit m T, where m is the ψ mass, leading in Fourier space to the following retarded self-energy of the mediator: Dropping the vacuum part and integrating over k 0 one obtains: The result so far is exact up to the fact that we neglected the ψ mass and ignored the vacuum contribution. Now, the Hard-Thermal-Loop approximation [66] assumes that the external energy p 0 and momentum |p| are smaller compared to the typical loop momentum |k| which is of the order temperature (hard), since the integrand contains n F . Expanding the term in the brackets to leading order in p 0 /|k| and |p|/|k|, all remaining integrals can be performed analytically leading to the finite result: where the Debye mass is defined as m 2 D = g 2 ψ T 2 /3. This result coincides with the result obtained from the imaginary time formalism, where instead one has to perform a sum over Matsubara frequencies.

Appendix D: Vacuum limit of two-particle spectral function
Here we will derive Eq. (VI.12) starting from Eq. (V.44) in the limit of vanishing finite temperature corrections. As we have briefly mentioned in Section V D, for the vacuum limit one has to carefully take into account the imaginary part i in the retarded equation, representing the small width which will be taken to be zero in the end. We will see that the result does not depend on this i -prescription as long as is small enough.
Suppose that the potential almost vanishes for a large enough r. For a Yukawa type potential, this is true for m V r 1. The parameter should be much smaller than this mass parameter, namely m V . In the case of the Coulomb potential one may introduce another small mass parameter to the gauge boson. In the end of the computation one can take it to be zero while keeping m V . For m V r 1, the homogeneous solutions, g >,< , can be well approximated by the plane wave: The Wronskian tells us that there exists a non-trivial relation between two coefficients C >,< . Since the Wronskian, W(r) = g > (r)g < (r) − g > (r)g < (r), does not depend on r, one may equate it at r = 0 and r → ∞, which yields We have taken to be zero in the end of the computation. Let us evaluate the integral given in Eq. (V.44) by means of Eq. (D.2). The first observation is that there is no imaginary part for g < if is zero from the beginning. Thus, the integrand becomes relevant only after r where we have used E = Mv 2 rel /4. Now we are in a position to discuss its relation to the conventional definition of the Sommerfeld enhancement factor. In the limit of → 0, the wave function propagates to infinity. Then one may obtain the Sommerfeld enhancement factor by extracting the amplitude of the wave function at the infinity, which is nothing but the relation S(v rel ) = |C > | 2 , see Ref. [41] for instance. Utilizing this relation, we finally get Eq. (VI.12). For conventional reason, let us give the s-wave Sommerfeld enhancement factor for the Coulomb case consistent with our equations: For the bound state, it is much easier to solve the equation directly rather starting from Eq. (V.44). One may express the spectral function by means of the wave functions for the bound states [34]: where ψ (B) n represents the normalized wave function for the n-th bound state. For instance, the wave function for the lowest energy state n = 1 is given by The decay rate of the bound state is related to its wave function at the origin. For the lowest state, one can easily show this from Eq. (D.9): (σv rel ) ψ (B) 1 (0) where the decay rate of the lowest bound state is given by (D.11) Similar calculation holds in the limit of negligible thermal width but finite real-part corrections. For this case one should substitute the kinetic energy E → E − V eff (∞) = Mv 2 rel /4 and similar for the bound state energy. Also one has to take smaller than m D and m V . thermodynamic picture like the grand canonical ensemble. It requires a rigorous proof of this claim, which we would like to give somewhere else.
Instead, we would like to give some approximations in order to obtain an analytic solution of Eq. (E.1). We restrict the discussion by assuming the system is in a grand canonical state. Utilizing the KMS condition for finite chemical potential, one can formally solve G +− η in terms of spectral function: In the last equality we assumed the DM gas to be dilute and approximated the Fermi-Dirac distribution as a Maxwellian. By KMS relation we have formally solved for the number density in terms of chemical potential and spectral function. The spectral function can be computed from the retarded component G R η , according to G ρ η,0 (p, ω) = G R η,0 (p, ω) − G A η,0 (p, ω) (see Section II). The EoM of the retarded correlator can be obtained from Eq. (E.1) by subtracting +− from the ++ component, e.g. G R η = G ++ η − G +− η . In subsequent sections we solve the retarded equation in various approximations, compute the spectral function and finally evaluate Eq. (E.3).

Ideal gas approximation
The ideal gas approximation can be defined as the zeroth order contribution in Eq. (E.1). This means we have to know the solution of the free retarded correlator. It can be obtained from the differential form of the retarded equations. In Fourier space of the microscopic Wigner coordinates it is given by G R η,0 (p, ω) = iδ i j /(ω − p 2 /(2M) + i ), leading to G ρ η,0 (p, ω) = δ i j (2π)δ(ω − p 2 /2M). Inserting this into Eq. (E.3) results in the ideal gas approximation of the free number density: We can invert this relation to obtain the chemical potential as a function of the number density in the ideal gas approximation: And similar expressions for anti-particle ξ.