Dark Matter Sommerfeld-enhanced annihilationand Bound-state decay at finite temperature

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

Tobias Binder Institute for Theoretical Physics, Georg-August University Göttingen, Friedrich-Hund-Platz 1, Göttingen, D-37077 Germany    Laura Covi Institute for Theoretical Physics, Georg-August University Göttingen, Friedrich-Hund-Platz 1, Göttingen, D-37077 Germany    Kyohei Mukaida Deutsches Elektronen-Synchrotron (DESY), Notkestraße 85, Hamburg, D-22607 Germany
July 27, 2019
###### Abstract

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.

preprint: DESY 18-123

## 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 Ade et al. (2016); Aghanim et al. (2018), 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 , set the most stringent bounds on the self-scattering cross section to be less than in the bullet cluster Randall et al. (2008) (in order to guarantee the production of elliptical halos Miralda-Escude (2002); Peter et al. (2013)). On the other hand, a DM self-scattering cross section of the order on dwarf-galactic scales, where velocities are of the order , 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 Tulin and Yu (2018).

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 Lee and Weinberg (1977); Ellis et al. (1984); Arcadi et al. (2018), such that the relic abundance in the Early Universe is set by the thermal 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 Ade et al. (2016); Aghanim et al. (2018). If the light mediators induce an attractive force between the annihilating DM particles, the total cross section is typically enhanced Hisano et al. (2004, 2007) which is often referred as the Sommerfeld enhancement Sommerfeld () or Sakharov enhancement Sakharov (1948). Additionally, such attractive forces can lead to the existence of DM bound-state solutions Pospelov and Ritz (2009); March-Russell and West (2009); Shepherd et al. (2009). This opens the possibility for conversion processes between scattering and bound states via radiative processes, influencing the evolution of the abundance of the stable scattering states during the DM thermal history. DM scenarios with Sommerfeld enhancement or with bound states have been widely studied in the literature Belotsky et al. (2008); Arkani-Hamed et al. (2009); Slatyer (2010); Feng et al. (2010); von Harling and Petraki (2014); Petraki et al. (2015); An et al. (2016); Asadi et al. (2017); Johnson et al. (2016); Hryczuk et al. (2011); Beneke et al. (2016); Freitas (2007); Hryczuk (2011); Harigaya et al. (2014); Harz et al. (2015); Beneke et al. (2015); Ellis et al. (2016); Kim and Laine (2017); El Hedri et al. (2017); Liew and Luo (2017); Mitridate et al. (2017); Harz and Petraki (2018a); Petraki et al. (2017); Harz and Petraki (2018b); Cirelli et al. (2007); March-Russell et al. (2008); Hisano et al. (2005); Cirelli et al. (2008, 2017); Fan and Reece (2013); Bhattacherjee et al. (2014); Ibe et al. (2015); Beneke et al. (2017); Berger et al. (2008); Blum et al. (2016); Biondini and Laine (2017, 2018); Biondini (2018) and it has been shown that the main effect of such corrections is to shift in the parameter space the upper bounds on the DM mass, otherwise the theoretically predicted DM density would be too large (overclosure bound).

Classic WIMP candidates with large corrections via Sommerfeld enhancement or bound states are particles charged under the electroweak interactions, like the Wino neutralino in supersymmetric models Jungman et al. (1996) or the first Kaluza-Klein excitation of the gauge boson in models with extra dimensions Servant and Tait (2003). For the supersymmetric case it was realized very early on by Hisano et al. (2004, 2007) that the Sommerfeld effect reduces the Wino density up to 30 % and pushes the mass of Wino Dark Matter candidate to few TeVs in order to obtain the correct relic density. These studies have later been extended to the case of general components of neutralino Hryczuk et al. (2011); Beneke et al. (2016). Similar and even stronger effects from the Sommerfeld enhancement and bound states were found in the case of coannihilation of the WIMP particle with charged or colored states Freitas (2007); Hryczuk (2011); Harigaya et al. (2014); Harz et al. (2015); Beneke et al. (2015); Ellis et al. (2016); Kim and Laine (2017); El Hedri et al. (2017); Liew and Luo (2017); Mitridate et al. (2017); Harz and Petraki (2018a). If the electroweakly charged Dark Matter is sufficiently heavy, the Sommerfeld enhancement or the presence of bound states due to the exchange of electroweak gauge or Higgs bosons, see e.g. Petraki et al. (2017); Harz and Petraki (2018b), are very generic as it was shown for example in the minimal Dark Matter model Cirelli et al. (2007); Mitridate et al. (2017) and in Higgs portal models March-Russell et al. (2008). In these cases, long-range force effects play an important role also for the indirect detection limits Hisano et al. (2005); Cirelli et al. (2008); Pospelov and Ritz (2009); March-Russell and West (2009); Cirelli et al. (2017) and especially for the Wino the Sommerfeld enhancement has lead to the exclusion of most parameter space Fan and Reece (2013); Bhattacherjee et al. (2014); Ibe et al. (2015); Beneke et al. (2017). Note that this effect can be important also when the Dark Matter is not itself a WIMP, but it is produced by WIMP decay out-of-equilibrium, like in the SuperWIMP mechanism Covi et al. (1999); Feng et al. (2003). Indeed in such scenario, the DM inherits part of the energy density of the mother particle and so any change in the latter freeze-out density is directly transferred to the superweakly interacting DM and can relax the BBN constraints on the mother particle Berger et al. (2008).

While a lot of effort has been made to compute quantitatively the effects of a 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 Burnier et al. (2008). For DM or heavy quark systems, the Sommerfeld enhancement at finite temperature has been discussed in the literature, where the chemical equilibration rate is i) estimated from linear response theory Bodeker and Laine (2013); Kim and Laine (2016, 2017) and ii) based on classical rate arguments Bodeker and Laine (2012), is then inserted into the 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 Biondini and Laine (2017, 2018); Biondini (2018). However, strictly speaking, the linear response theory method is only applicable if the system is close to chemical equilibrium. Indeed the computation has been done in the 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.

## Ii Real-time formalism prerequisites

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 prescription111In ordinary QFT the initial vacuum state appearing in correlation functions is equivalent up to a phase to the final vacuum state. For this special situation the operators are ordered along the ’flat’ time axis ranging from to . By means of LSZ reduction formula it is then possible to relate correlation functions to the matrix and compute cross sections. This 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 Stefanucci and van Leeuwen (2013). of the CTP contour, as illustrated in Fig. 1. The time contour in the Keldysh-Schwinger prescription consists of two branches denoted by and . The upper time contour ranges from the initial time to while the lower contour is considered to go from back to . Therefore, times on the branch are said to be always later compared to the times on . The time ordering of operator products on can generically be written as

 TC[O1(t1)...On(tn)]≡∑P(−1)F(p)θC(tp(1),tp(2))...θC(tp(n−1),tp(n))Op(1)(tp(1))...Op(n)(tp(n)), (II.1)

where the sum is over all the permutations of the set of operators and is the number of permutations of fermionic operators. The unit step function and the delta distribution on the Keldysh-Schwinger contour is defined as

 θC(t1,t2)≡⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩θ(t1−t2) if t1,t2∈τ+θ(t2−t1) if t1,t2∈τ−1 if t1∈τ−,t2∈τ+0 if t1∈τ+,t2∈τ−^=(θ(t1−t2)01θ(t2−t1)),δC(t1,t2)^=(δ(t1−t2)00−δ(t1−t2)). (II.2)

Correlation functions, i.e. contour -ordered operator products averaged over all states where the weight is the density matrix of the system denoted by , are defined by

 (II.3)

Let us introduce commonly used notations and properties of 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 and on and hence four different components of a general two-point function denoted by , where in matrix notation it can be written as

 G(x,y)≡⟨TCψ(x)ψ†(y)⟩=θC(x0,y0)⟨ψ(x)ψ†(y)⟩∓θC(y0,x0)⟨ψ†(y)ψ(x)⟩^=(G++(x,y)G+−(x,y)G−+(x,y)G−−(x,y)). (II.4)

Here, means and with for and the four different components of are defined as:

 G−+(x,y) ≡⟨ψ(x)ψ†(y)⟩, (II.5) G+−(x,y) ≡∓⟨ψ†(y)ψ(x)⟩, (II.6) G++(x,y) ≡θ(x0−y0)G−+(x,y)+θ(y0−x0)G+−(x,y), (II.7) G−−(x,y) ≡θ(x0−y0)G+−(x,y)+θ(y0−x0)G−+(x,y), (II.8)

where on the r.h.s. of the second line applies for fermionic (bosonic) field operators. From these definitions one can recognize that not all components are independent, namely the following relation holds

 G++(x,y)+G−−(x,y)=G+−(x,y)+G−+(x,y). (II.9)

Furthermore, let us introduce retarded and advanced correlators defined by

 GR(x,y) ≡θ(x0−y0)[G−+(x,y)−G+−(x,y)]=G++(x,y)−G+−(x,y)=−G−−(x,y)+G−+(x,y), (II.10) GA(x,y) ≡−θ(y0−x0)[G−+(x,y)−G+−(x,y)]=G++(x,y)−G−+(x,y)=−G−−(x,y)+G+−(x,y), (II.11)

as well as the spectral function given by:

 Gρ(x,y)≡GR(x,y)−GA(x,y)=G−+(x,y)−G+−(x,y). (II.12)

From these definitions we can derive further useful properties:

 GA(x,y) =−[GR(y,x)]†,G+−(x,y)=[G+−(y,x)]†,G−+(x,y)=[G−+(y,x)]†. (II.13)

In the case of free (unperturbed) propagators , the following semigroup property holds:

 GR0(x,y) =∫d3zGR0(x,z)GR0(z,y),for tx>tz>ty. (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 of products of correlators can be written as:

 C++(x,y) =[∫z0∈Cdz0∫d3zA(x,z)B(z,y)]x0,y0∈τ+=[(∫τ+% dz0+∫τ−dz0)∫d3zA(x,z)B(z,y)]x0,y0∈τ+ =∫∞−∞dz0∫d3zA++(x,z)B++(z,y)+∫−∞∞dz0∫d3zA+−(x,z)B−+(z,y) =∫∞−∞d4z(A++(x,z)B++(z,y)−A+−(x,z)B−+(z,y)). (II.15)

Eq. (II.15) is called a Lagereth rule and it is straightforward to work out similar rules for, e.g. different components or double integrations of 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

 T =(tx+ty)/2,t=tx−ty, (II.16) R =(x+y)/2,r=x−y. (II.17)

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(txx,tyy). (II.18)

In all computations, the tilde will be dropped such that we can write for the Fourier transformation of :

 G(ω,p,R,T)=∫dtd3rei(ωt−p⋅r)G(t,r,R,T). (II.19)

One of the great advantages of separating microscopic () and macroscopic () variables according to the Wigner transformation is that Fourier transformations of integral expressions can be considerably simplified by using the gradient expansion. For example, Eq. (II.15) in Fourier space can be written as

 C++(ω,p,R,T) (II.20) ≃A++(ω,p,R,T)B++(ω,p,R,T)−A+−(ω,p,R,T)B−+(ω,p,R,T), (II.21) GAB ≡e−i(∂AT∂Bω−∂Aω∂BT−∇AR⋅∇Bp+∇Ap⋅∇BR)/2≃1. (II.22)

Throughout this work, such exponentials containing derivatives are always approximated as in the last line, defining the leading order term of the gradient expansion. For homogeneous and isotropic systems the correlators do not depend on and thus for the spatial part the leading order term is exact. For a discussion of the validity of the leading order term of the temporal part we refer to Beneke et al. (2014). Let us emphasize here, that for typical DM scenarios the leading order term is always assumed to be a very good approximation.

Next, important properties of 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

 ^ρ∝e−βH, (II.23)

where is the Hamiltonian of the system and factor is the inverse temperature of the system. The density matrix in thermal equilibrium can be formally seen as a time evolution operator, where the inverse temperature is regarded as an evolution in the imaginary time direction. Making use of the cyclic property of the trace it can be shown that under the assumption of equilibrium the components are related via

 G−+(t1−t3)=∓G+−(t1−t3+iβ). (II.24)

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:

 G−+(ω,p)=∓eβωG+−(ω,p). (II.25)

From this condition, various equilibrium relations follow:

 G+−(ω,p) (II.26) G++(ω,p) =GR(ω,p)+GA(ω,p)2+[12∓nF/B(ω)]Gρ(ω,p). (II.27)

where the Fermi-Dirac or Bose-Einstein phase-space densities are given by . Thus in equilibrium, all correlator components can be calculated from the retarded/advanced components, where the spectral function is related to those via Eq. (II.12).

General 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 . 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 at the initial time that we take to in the Schwinger-Keldysh prescription of the contour, and the first one takes the average of with the weight of . Now, to derive the EoM for two-point correlators, let us consider an infinitesimal perturbation , satisfying . By relying on the measure-invariance principle under this transformation, one obtains the EoM of the two-point correlators from

 ⟨O′†(y)⟩ =⟨O†(y)⟩+ϵ(x)+i∫x∈Cϵ(x)⟨TCδSδLO†(x)O†(y)⟩+O(ϵ2) (II.29) ⇒0 =δC(x,y)+i⟨TCδSδLO†(x)O†(y)⟩, (II.30)

where represents a derivative acting from the left. The same procedure can be applied to the case of , as well as for deriving the EoM of higher correlation functions. The relation between the abstract EoM of correlators and differential equations for observables will be part of the next Section. In general, a correlator EoM depends on higher and lower correlators which is called the 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 Beneke et al. (2014). 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 two- and four-point 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:

 L=¯χ(i⧸∂−M)χ+gχ¯χγμχAμ+¯ψ(i⧸∂−m)ψ+gψ¯ψγμψAμ−14FμνFμν. (III.1)

The particles of the equilibrated plasma environment with temperature are the abelian mediators and the light fermionic particles with mass . Fermionic DM is assumed to be nonrelativistic, i.e. . All fermionic particles are considered to be of Dirac type. We assume the mediator to be massless, however, we provide the final results also for the case of a massive with mass .

Let us illustrate how we can get the DM nonrelativistic effective action in the thermal medium of light particles. It is obtained in two steps. First, hard modes of are integrated out. In this limit, the DM four component spinor splits into two parts, a term for the particle denoted by the 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:

 SNR[η,ξ]=∫\displaylimitsx∈Cη†(x)[i∂t+∇22M]η(x)+ξ†(x)[i∂t−∇22M]ξ(x)+∫\displaylimitsx,y∈Cig2χ2J(x)D(x,y)J(y)+iO†s(x)Γs(x,y)Os(y). (III.2)

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

 J(x)≡η†(x)η(x)+ξ†(x)ξ(x),D(x,y)≡⟨TCA0(x)A0(y)⟩, (III.3)

respectively. The last term in Eq. (III.2) describes the annihilation part and we only keep the -wave contribution

 Os(x)≡ξ†(x)η(x),Γs(x,y)≡π(α2χ+αχαψ)M2(δ4(x−y)02δ4(x−y)δ4(x−y)), (III.4)

with the fine structure constant being , and summation over the spin indices are implicit. is shown in the matrix representation of the CTP formalism, see e.g. Eq. (II.4) in previous Section. Hence the delta functions on the right-hand-side are defined on the usual real-time axis. Similar to the vacuum theory, the annihilation part can be computed by cutting the box diagram (containing two ) and the vacuum polarization diagram (containing one and a loop of light fermions ), where now all propagators are defined on the CPT-contour. Finite temperature corrections to these hard processes in are neglected222Assuming free plasma field correlators in the computation of is a good approximation since the energy scale of the hard process is which is much larger for nonrelativistic particles than typical finite temperature corrections being of the order . Consequently, the dominant thermal corrections should be in the modification of the long-range force correlator , where the typical DM momentum-exchange scale enters which is much lower compared to the annihilation scale. and for a derivation we refer to Appendix B.

In our effective action Eq. (III.2), we have discarded higher order terms in (like magnetic interactions) and also interactions with ultra-soft gauge bosons,333 To fully study the dynamics of the bound state formation Brambilla et al. (2008); von Harling and Petraki (2014); Petraki et al. (2015, 2017) at late times of the freeze-out process (e.g. with 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.4. Note that at high enough temperature those processes are typically efficient, leading just to the ionization equilibrium among bounded and scattering states. To estimate the time when the ionization equilibrium is violated concretely, we have to take into account these processes in the thermal plasma, which will be presented elsewhere. since we focus on threshold singularities of annihilations at the leading order in the coupling and the velocity . Furthermore, our effective field theory is 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 Hisano et al. (2004, 2005, 2007) and belongs to the ++ component of . Thus, as a first result we have generalized the annihilation term towards the CTP contour. Another one stems from the gauge boson propagator that encodes interactions with the thermal plasma. While the annihilation term containing in our action breaks the number conservation of DM, the interaction term containing can not. From this observation, one may anticipate that the 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 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 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.1. Second, in Section III.2 the kinetic equations for the DM correlators are derived. We show how to extract from these equations the number density equation, including finite temperature corrected processes for the negative energy spectrum (bound-state decays) as well as for the positive energy contribution (Sommerfeld-enhanced annihilation) in one single equation.

### iii.1 Thermal corrections to potential term

In this section, we briefly summarize how the electric component of the mediator propagator gets modified by the thermal presence of ultrarelativistic fields. The plasma environment is regarded to be perturbative and in the one-particle self-energy framework we can write down the Dyson equation on the Schwinger-Keldysh contour for the mediator:

 Dμν(x,y) =D0μν(x,y)−i∫w,z∈CD0μα(x,w)Παβ(w,z)Dβν(z,y), (III.5) Πμν(x,y) =(−i)g2ψ(−1)Tr[γμS(x,y)γνS(y,x)], (III.6)

where are unperturbed correlators. and are assumed to be in equilibrium and thus, according to the discussion in Section II, we only need to compute the retarded/advanced propagators. From those, we can construct all other components by using KMS condition. The Dyson equation for retarded (advanced) 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:

 ΠμνR(P)=Πμν++−Πμν+−=ig2ψ∫d4K(2π)4(Tr[γμS++(K−P)γνS++(K)]−Tr[γμS+−(K−P)γνS−+(K)]). (III.8)

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 Braaten and Pisarski (1990) to extract leading thermal corrections.444 Let us briefly summarize here the assumptions of the HTL approximation. First, we drop all vacuum contributions and only keep temperature dependent parts. Second, we assume the external energy and momentum to be much smaller than typical loop momentum which is of the order temperature (see Appendix C). The discussion of the validity of the HTL approximation depends on where the dressed mediator correlator is attached to. One can not naively argue for the case where one would attach it to the DM correlator that the external momentum is the DM momentum which is of course much larger than temperature, thus invalidating the HTL approximation by this argumentation. For example, in our case the dressed mediator correlator enters the DM single-particle self-energy (see Fig. 2 and Section V.3) and the dominant energy and momentum region in the loop diagram is where HTL effective theory is valid. In the HTL approximation we are allowed to resum the 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 of the mediator propagator in the HTL resummed self-energy approximation and in the Coulomb gauge the following:

 DR,A(ω,p) =−iω2−p2+Π00R,A(ω,p)±isign(ω)ϵ, (III.9) Π00R,A(ω,p) =−m2D[1−ω2|p|ln(∣∣∣ω+|p|ω−|p|∣∣∣)±iω|p|π2θ(|p|2−ω2)], (III.10)

where we introduced the Debye screening mass . One can recognize that there is correction to the real part of the mediator propagator as well as a branch cut for space-like exchange. Using the equilibrium relation Eq. (II.27), the component of in the static limit reads

 limω→0D++(ω,p)=limω→0[DR(ω,p)+DA(ω,p)2+(12+nB(ω))Dρ(ω,p)]=ip2+m2D+πT|p|m2D(p2+m2D)2, (III.11)

while for a massive mediator we have simply

 limω→0D++(ω,p)=ip2+m2V+m2D+πT|p|m2D(p2+m2V+m2D)2. (III.12)

The static 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 Bellac (2011). 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 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.

### iii.2 Exact DM number density equation from correlator equation of motion

The main purpose of this section is to derive the equation for the DM number density directly from the exact EoM of our nonrelativistic action. Defining the DM particle and anti-particle correlators as

 Gη(x,y) ≡⟨TCη(x)η†(y)⟩, (III.13) Gξ(x,y) ≡⟨TCξ(x)ξ†(y)⟩, (III.14)

we derive respective EoM from the path-integral formalism, as briefly explained at the end of Section II, for the nonrelativistic effective action given in Eq. (III.2):

 [i∂x0+∇2x2M]Gη(x,y) =iδC(x,y)−ig2χ∫\displaylimitsz∈CD(x,z)⟨TCη(x)J(z)η†(y)⟩−i∫\displaylimitsz∈CΓs(x,z)⟨TCξ(x){ξ†(z)η(z)}η†(y)⟩, (III.15) [−i∂y0+∇2y2M]Gη(x,y) =iδC(x,y)−ig2χ∫\displaylimitsz∈CD(y,z)⟨TCη(x)J(z)η†(y)⟩−i∫\displaylimitsz∈CΓs(z,y)⟨TCη(x){η†(z)ξ(z)}ξ†(y)⟩, (III.16) [i∂x0−∇2x2M]Gξ(x,y) =iδC(x,y)−ig2χ∫\displaylimitsz∈CD(x,z)⟨TCξ(x)J(z)ξ†(y)⟩−i∫\displaylimitsz∈CΓs(x,z)⟨TCη(x){η†(z)ξ(z)}ξ†(y)⟩, (III.17) [−i∂y0−∇2y2M]Gξ(x,y) =iδC(x,y)−ig2χ∫\displaylimitsz∈CD(y,z)⟨TCξ(x)J(z)ξ†(y)⟩−i∫\displaylimitsz∈CΓs(y,z)⟨TCξ(x){ξ†(z)η(z)}η†(y)⟩. (III.18)

The anticipated structure in Eq. (III.15)-(III.18) shows the dependence of the two-point correlators on higher correlation functions. Here, the curly brackets stand for the summation over the spin indices and is the current as already introduced in Eq. (III.3). It might be helpful to mention that we used a special property of two-point functions of Hermitian bosonic field operators: . This exact property can be verified directly from the definition in Eq. (III.3).

In the following, the number density equation of particle and 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 in Eq. (III.2). For this purpose, let us switch off the annihilation term in and seek for conserved quantities. In this limit, the theory has the following global symmetries: and . The associated Noether currents for DM particle and anti-particle are:

 Jμη(x)=(η†(x)η(x),12iMη†(x)↔∇η(x)),Jμξ(x)=(ξ(x)ξ†(x),12iMξ(x)↔∇ξ†(x)). (III.19)

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 . For the particle DM, we obtain as an intermediate result after all these steps:

 i∂μ⟨Jμη⟩= −(i∂x0+i∂y0+∇2x2M−∇2y2M)TrGη(x,y)∣∣y→x = +ig2χ∫z∈CD(x,z)[⟨TC{η(x)η†(x)}J(z)⟩−⟨TC{η(x)η†(x)}J(z)⟩] (III.20)

The trace and the curly brackets indicate the summation over the spin indices. It is important to note that the first line in the second equality cancels out, even in the case of a fully interacting correlator including finite temperature corrections. Thus, we confirm from the EoM that, e.g. non-Hermitian potential corrections arising from soft thermal DM- scatterings in the HTL approximation of [see Eq. (III.11)], never violate the current conservation of each individual DM species. For a homogeneous and isotropic system (vanishing divergence of current density) this would mean that the individual number densities of particles and antiparticles do not change by self-scattering processes, real physical DM- scatterings, soft DM- scatterings or other finite temperature corrections leading to potential contributions in .

It can be recognized, that the current conservation is only violated by the annihilation term in the last line in Eq. (III.20), since this contribution does not cancel to zero. This term can be simplified by using Eq. (III.4) and by fixing the time component to either or . We have explicitly checked that both choices of lead to the same final result. With the definition of the four-point correlator on the closed-time-path contour

 Gηξ,s(x,y,z,w)≡⟨TCηi(x)ξ†i(y)ξj(w)η†j(z)⟩, (III.21)

we obtain our final form of the current equations.

Current equations for particle and anti-particle :

 ∂μ⟨Jμη(x)⟩ =−2π(α2χ+αχαψ)M2G++−−ηξ,s(x,x,x,x), (III.22) ∂μ⟨Jμξ(x)⟩ =−2π(α2χ+αχαψ)M2G++−−ηξ,s(x,x,x,x). (III.23)

The current conservation is only violated by contributions coming from . This is consistent with the expectations from the symmetry properties of the action when annihilation is turned on. Namely, only a linear combination of both global transformations leaves the action invariant which leads to the conservation of , which is nothing but the DM asymmetry current conservation. The conservation of the total DM number density, , is violated by the annihilation term.

Before we discuss the 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 . 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 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:

 ˙nη+3Hnη =−2(σvrel)[G++−−ηξ,s(x,x,x,x)−G++−−ηξ,s(x,x,x,x)∣∣eq], (III.24) ˙nξ+3Hnξ =−2(σvrel)[G++−−ηξ,s(x,x,x,x)−G++−−ηξ,s(x,x,x,x)∣∣eq]. (III.25)

The tree-level s-wave annihilation cross section of our system was defined as and in the last term means the evaluation at thermal equilibrium. Note that in the CTP formalism a cross section strictly speaking does not exist. The reason why this result is equal to the vacuum computation is because we computed the annihilation part at the leading order, where it is expected that zero and finite temperature results should coincide. The correlation function however is fully interacting. We summarize with two concluding remarks on our main results:

• 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 . 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 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, contains both, the bound state and scattering state contributions at the same time and at finite temperature they turn out to be not separable as it is sometimes done in vacuum computations. Bound state contributions will automatically change the cross section into a decay width and thus, appearing on the l.h.s. of Eq. (III.24) is the total number of particles including the ones in the bound states and similar interpretation for the anti-particle .

• Particle number-conservation: In Section III.1, we have seen that the thermal corrections to the mediator propagator can contain, next to the real Debye mass, an imaginary contribution. It was shown that these 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 , the nonrelativistic action in Eq. (III.2) has two global symmetries and . The conserved quantities are the particle and anti-particle currents in Eq. (III.22) and (III.23) in the limit (vanishing r.h.s). When annihilation is included, the nonrelativistic action is only invariant if both global transformations are performed at the same time, resulting in the conserved asymmetry current . We conclude that thermal corrections can never violate these symmetries, even not at higher loop level. On the other hand, the solution of the Sommerfeld factor is contained in and hence the annihilation rate will depend on the thermal loop corrected longe-range mediator , 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 . 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 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 are given by:

 [i∂x0+∇2x2M]Gη(x,y) =iδC(x,y)−ig2χ∫z∈CD(x,z)[Gηξ(x,z,y,z)−Gηη(x,z,y,z)], (IV.1) [i∂x0+∇2x2M]¯Gξ(x,y) =iδC(x,y)−ig2χ∫z∈CD(x,z)[Gηξ(z,x,z,y)−Gξξ(x,z,y,z)], (IV.2) [i∂x0+∇2x2M]Gηξ(x,y,z,w) =iδC(x,z)¯Gξ(y,w)−ig2χ∫¯x∈CD(x,¯x)⟨TCη(x)J(¯x)ξ†(y)ξ(w)η†(z)⟩, (IV.3)

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

 ¯Gξ(x,y) ≡⟨TCξ†(x)ξ(y)⟩, (IV.4) Gηξ(x,y,z,w) ≡⟨TCη(x)ξ†(y)ξ(w)η†(z)⟩, (IV.5) Gηη(x,y,z,w) ≡⟨TCη(x)η(y)η†(w)η†(z)⟩, (IV.6) Gξξ(x,y,z,w) ≡⟨TCξ†(x)ξ†(y)ξ(w)ξ(z)⟩. (IV.7)

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 , 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 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δC(y,w)[Gηξ(x,¯x,z,¯x)−Gηη(x,¯x,z,¯x)]−iδC(¯x,y)Gηξ(x,y,z,w), (IV.8)

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