Re-derived overclosure bound for the inert doublet model

We apply a formalism accounting for thermal effects (such as modified Sommerfeld effect; Salpeter correction; decohering scatterings; dissociation of bound states), to one of the simplest WIMP-like dark matter models, associated with an"inert"Higgs doublet. A broad temperature range T ~ M/20...M/10^4 is considered, stressing the importance and less-understood nature of late annihilation stages. Even though only weak interactions play a role, we find that resummed real and virtual corrections increase the tree-level overclosure bound by 1...18%, depending on quartic couplings and mass splittings.


Introduction
Tight constraints from the LHC and from direct and indirect detection experiments have put many simple dark matter models under tension in recent years. This calls for new ideas in model building, but perhaps also for new precision in the computations on which a given dark matter scenario is based. Indeed, as the LHC pushes up the dark matter mass scale, it also increases the temperature at which dark matter density was fixed. Then, however, Standard Model weak interactions, which play a role in most dark matter computations, can be modified by thermal effects. If the freeze-out temperature is T > ∼ 160 GeV, the Higgs mechanism "melts away" [1], whereby weak interactions display phenomena normally only associated with strong interactions.
The purpose of this paper is to present a step-by-step implementation of a formalism which can account for relevant thermal effects [2], 1 and whose principal applicability has been tested against non-perturbative lattice simulations by using the annihilation of heavy quarks in QCD as an analogue for dark matter annihilation [9]. Among our goals are to check whether thermal modifications affect the well-known Sommerfeld enhancement (cf. e.g. refs. [10][11][12][13]), and how to include the classic Salpeter correction (cf. e.g. ref. [14]).
The premise of the framework is to make use of a heavy-mass or "non-relativistic" expansion for the dark matter particles. Given that in the classic WIMP paradigm dark matter gradually freezes out at a temperature T ∼ M/20...M/10 4 , where M is the dark matter mass scale, there should be no doubt about the validity of this approximation.
Within the non-relativistic regime, the framework accounts for a number of thermal effects, such as that the vacuum masses of W ± , Z 0 are replaced by thermal Debye masses as the temperature increases; that the weak mixing angle evolves with the temperature; that weak interactions mediate fast scatterings of the dark matter particles, transforming them into each other and thereby affecting the nature of their annihilation process; that similar interactions also change the effective mass of the dark matter particles through the Salpeter correction; and that in some cases dark matter particles can form bound states. As far as the coannihilation of non-degenerate dark matter particles goes, the formalism can also be nicely contrasted with the classic Boltzmann equation approach of ref. [15].
To put the study in context, we remark that there has been recent interest in including next-to-leading order (NLO) corrections into dark matter computations. Here we are more concerned with the fact that most computations are formally incomplete even at leading order (LO), as far as near-threshold thermal effects go [2]. In principle, the inclusion of NLO corrections is also possible within the same formalism, notably by adding operators suppressed by ∼ ∇ 2 /M 2 to eq. (2.3) and NLO corrections to the coefficients given in eqs. (2.4)-(2.6), however this is not pursued here.
The model with which we choose to illustrate the formalism is a simple extension of the Standard Model through an additional "inert" Higgs doublet [16][17][18]. Many dark matter computations have been carried out for various parameter corners of this model (cf. e.g. refs.  and references therein; we particularly recommend ref. [25] for a general overview), and our conclusions do not differ qualitatively from these, even though visible effects from hitherto unconsidered processes can be observed.
The plan of this paper is the following. After introducing the 4-particle operators that mediate dark matter decays in the heavy-mass limit (sec. 2), we recall how they determine the thermal dark matter annihilation rate (sec. 3). Subsequently the key tools of the formalism, namely time-dependent medium-modified Schrödinger equations governing the "slow" dynamics within the dark sector, are elaborated upon (sec. 4). After presenting numerical solutions and the overclosure bound (sec. 5), we turn to conclusions and an outlook (sec. 6).

4-particle operators
In the inert doublet model (IDM), the Standard Model is supplemented by an additional Higgs doublet, χ, which does not couple to fermions because of an unbroken discrete Z(2) symmetry. Denoting by φ the Standard Model Higgs doublet and by D µ the corresponding covariant derivative, the Standard Model Lagrangian is modified by the additional terms (2.1) The notation λ 1 is reserved for the Standard Model Higgs self-coupling, If the mass scale M is much larger than the electroweak scale, M ≫ m W , the χ particles annihilate efficiently into the Standard Model ones. The annihilations that we are interested in happen in the temperature range T ∼ M/20...M/10 4 , in which the average velocity is v ∼ T /M ≪ 1. Therefore the annihilating particles are non-relativistic. Non-relativistic annihilations can be described by 4-particle operators, arranged as an expansion in 1/M 2 [42]. If we write non-relativistic on-shell fields in terms of annihilation and creation operators as 2) then at leading order in 1/M 2 there are four "absorptive" operators that play a role: 2 Here sums over the isospin components p, q, r, s ∈ {1, 2} are implied, and T a ≡ σ a /2, where σ a are the Pauli matrices.
We have computed the coefficients c 1 , ..., c 4 in eq. (2.3) in general R ξ gauges at leading non-trivial order, verifying their gauge independence for ξ < M 2 /m 2 Z : 3 Here g 1 and g 2 are the U Y (1) and SU L (2) gauge couplings, respectively. The couplings should be evaluated at a renormalization scale ∼ 2M . The same values of the coefficients can be extracted from ref. [25].
If λ 4 = 0 or λ 5 = 0, or if Standard Model radiative corrections are considered, different components of χ are non-degenerate in mass. In this case the doublets C and D can be written as The operators in eq. (2.3) split into a 10 × 10 matrix in the field space of eq. (2.7), which is given (with a slightly different notation) in eqs. (4.25)-(4.28) below. 2 As is characteristic of an effective theory approach, there are in principle infinitely many higher-dimensional operators, suppressed by increasing powers of 1/M 2 . The four operators here are the only ones at order 1/M 2 . The coefficients of these operators contain both a real part and an imaginary (i.e. absorptive) part [42]. Only the imaginary parts are relevant for us [2]: in accordance with the optical theorem, they represent matrix elements squared of real processes in which the heavy particles annihilate into Standard Model ones. The annihilations are two-particle annihilations; therefore the matrix elements squared contain four field operators, two annihilation operators for a process, and two creation operators for its conjugate. In eqs. (2.4)-(2.6) the coefficients of these operators are given at leading order, corresponding to a tree-level annihilation cross section. One strength of the effective theory approach is that if needed, it would be fairly straightforward to compute NLO corrections to the coefficients. Even more importantly, soft thermal corrections to the annihilation processes (cf. fig. 1 for an illustration) can be included beyond a quasi-particle approximation, and up to the non-perturbative level in the case of strong interactions [9]. 3 For ξ ≫ M 2 /m 2 Z ≫ 1 the results change qualitatively and therefore unitary gauge is not viable.

Rate equations and effective cross sections
As discussed in ref. [15], the only physically reasonable "slow variable" of the problem at hand is the total number density of dark matter particles, n ≡ i=±,0,0 n i . Within a Boltzmann approach, ref. [15] established that n evolves according to the Lee-Weinberg equation [43,44], whereṅ is the covariant time derivative in an expanding background, and is an effective cross section for 2 → 2 annihilations from the dark sector. In our case the total equilibrium number density reads, at tree-level, and K 2 is a modified Bessel function. We note in passing that radiative corrections to eq. (3.3) can be determined as explained in ref. [2]. The most important is the so-called Salpeter correction, which modifies the rest mass of a non-relativistic particle by an amount ∆M T ∼ −α 3/2 T < 0, where α is a weak fine-structure constant (cf. e.g. ref. [5]). This is specified in more detail in sec. 5 (cf. eq. (5.6)).
In contrast to eq. (3.1), the formalism of ref. [2] takes as a starting point an equation based on general linear response theory, having thus the form [45] n = −Γ chem n − n eq + O n − n eq 2 , where Γ chem can be called the chemical equilibration rate. In the remainder of this paper, we wish to make close contact with standard literature, and therefore prefer to use the form of eq. (3.1). Linearizing eq. (3.1) in deviations from equilibrium leads us to identify In the absence of a first-principles argument beyond the linear response level, we rely on the form of eq. (3.1) on how first and higher order deviations are related to each other. The strength of the linear response approach is that it permits to relate the equilibration rate Γ chem to a correlator evaluated in equilibrium, without assuming weak interactions or the validity of a quasi-particle description necessary for a Boltzmann treatment [45]. Specifically, when the reactions responsible for equilibration are described by operators of the type in eq. (2.3), Γ chem is to first order proportional to the thermal expectation value of δL abs [9]. Inserting the proportionality coefficient and expressing the result through eq. (3.5), we obtain Because the annihilation operators are positioned to the right in eq. (2.3), the vacuum state does not contribute to the expectation value in eq. (3.6). Therefore γ i is exponentially suppressed by ∼ e −2M/T , with the Boltzmann factor cancelling against that from n 2 eq . Eq. (3.6) represents a generalization of eq. (3.2). The matrix structure of σ ij corresponds to matrix-like Schrödinger equations satisfied by the wave functions of the annihilating pair (cf. table 1), and the weights n eq i in eq. (3.2) correspond to threshold locations in the Laplace transform in eq. (4.39). At the same time eq. (3.6) goes beyond eq. (3.2) in several respects, for instance by permitting for a systematic inclusion of virtual thermal effects in the computation of individual cross sections, and also of real thermal scatterings of the dark matter particles off Standard Model particles, as discussed in more detail in secs. 4.1 and 4.2.

General goal and physical interpretation
In the notation of ref. [15], the cross sections in eq. (3.2) describe the processes where X, X ′ are Standard Model particles. These are "slow" processes: the likelihood that a dark matter particle finds a partner with which to annihilate is suppressed by a Boltzmann factor, so that the rate is Γ ∼ α 2 M 2 k e −E/T . However, the χ i -particles also experience "fast" reactions which have no Boltzmann suppression associated with them. These are of the type given in eqs. (6b) and (6c) of ref. [15]: These reactions keep the dark matter particles in kinetic equilibrium, and also change them into each other, guaranteeing chemical equilibrium within the dark sector, with each species contributing with its proper number density n i eq into eq. (3.2). If there are bound states in the dark sector, further "fast" processes can be added, notably where we assume that the binding energy is small, ∆E ∼ α 2 M < ∼ πT . Of course the same reactions are also present without bound states, and can change the annihilating pair into a different gauge or spin state. In addition, processes with virtual X exchange are important, leading e.g. to the Sommerfeld effect. The description based on eq. (3.6) goes beyond eq. (3.2) in that the indirect effect of the reactions in eqs. (4.2)-(4.5) can be included in a more "differential" form. Specifically, the fast reactions in eq. (4.2) give thermal masses to the dark matter particles, which change the kinematics of the reactions in eq. (4.1), leading e.g. to the Salpeter correction whereby the location of the 2-particle threshold gets modified. The fast reactions also induce thermal interaction rates, which decohere quantum-mechanical phases and thereby affect cross sections. Likewise the Sommerfeld effect and the possible emergence of bound states are included, through the solution of dynamical (time-dependent) Schrödinger equations. Thereby there is no need to assume the validity of a quasi-particle picture in the dark sector.

On the applicability of the Schrödinger description
Despite its strengths, an effective Schrödinger description as outlined in sec. 4.1 is only valid in a certain parametric regime. Indeed its justification requires an analysis of the different energy and momentum scales contributing to the problem. For near-threshold problems at finite temperature, several different scales play a role. A thermally modified Schrödinger approach in the form implemented below can be used for addressing energy scales ∆E ∼ α 2 M provided that (cf. e.g. refs. [46][47][48][49]) where α ∼ g 2 /(4π). In this situation the scale gT , which is the Debye scale representing typical energies/momenta of soft Standard Model excitations, can be integrated out, so that no Standard Model fields appear in the description of the "slow" dynamics. An example of an excitation associated with the scale gT is an electric dipole ∼ r · gE. As discussed in ref. [48], such dipoles cause transitions between pairs in different gauge representations, as appear in the operators of eq. (2.3). Specifically, integrating out the E fields and the pairs in repulsive channels generates a thermal interaction rate affecting the dynamics of the pair in an attractive channel [48]. Now, the interaction rate in the attractive channel is a slow rate: the annihilating pair is in a gauge-singlet state and only a dipole contribution is left over, Γ ∼ α 2 T 3 r 2 . Therefore, Γ can be part of an effective slow quantum-mechanical description.
In contrast, the generic interaction rates in the Standard Model, and in particular the interaction rates of the heavy χ pairs in gauge non-singlet channels, are of order αT . This is . . .  1), yielding a substantial below-threshold tail akin to that appearing below the top-antitop threshold in vacuum [50].
We have adopted a procedure here in which the contributions of the repulsive channels are estimated in two ways: either including the below-threshold tail, or omitting it. The difference of the results is used for estimating the theoretical uncertainties of our computation from thermal effects which are formally of NLO magnitude.
Having introduced the four-particle operators (cf. eq. (2.3)) and the Schrödinger approach, we can briefly comment on the different stages of the annihilation process. According to the scale hierarchy in eq. (4.6), there are two well-separated classes of processes: those occurring at the hard scale, M , and those typical of the soft scales, either thermal or non-relativistic (cf. fig. 1). The latter account for several interactions with particles from the heat bath which are resummed by a thermally modified Schrödinger equation. In the end the dark matter particles annihilate into Standard Model ones. This happens at a typical distance scale of order 1/M which is not resolved by the larger medium length scales. Hence an effective point-like interaction is responsible for the hard process. Such a factorization manifests itself in the effective cross section, eq. (5.3) below, where the hard coefficients from eqs. (2.4)-(2.6) multiply thermal expectation values capturing the soft physics.

Degenerate limit
We start by considering the degenerate limit, i.e. M ≡ M 0 = M0 = M ± . Each of the expectation values in eq. (3.6) can be expressed as a Laplace transform of a spectral function, denoted by ρ i (cf. eqs. (4.21)-(4.24) below). Under the assumptions discussed in sec. 4.2 and going over to non-relativistic center-of-mass coordinates, the spectral function is in turn an imaginary part of a Coulomb Green's function [2]: where N i is a normalization factor giving the number of contractions related to O i : In center-of-mass coordinates the Laplace transform reads where M ≫ Λ ≫ α 2 M is a cutoff restricting the average to the non-relativistic regime. According to eq. (3.6), the physical result is 4 i=1 c i γ i , with c i given in eqs. (2.4)-(2.6). In the free limit, V i → 0, the spectral function from eqs. (4.7) and (4.8) reads ρ Carrying out the Laplace transform in eq. (4.10), inserting n (0) eq = 4 M T 2π 3 2 e −M/T from eq. (3.3), and plugging into eq. (3.6), we obtain the value of σ eff v for a degenerate system and to leading order in 1/M 2 and α: In order to go beyond eq. (4.11), we include the potentials V i for the various channels in eq. (4.7). It is helpful to introduce the notation where ... T denotes a time-ordered propagator and the gauge potentials have been expressed with the sign conventions of the imaginary-time formalism. For instance (cf. appendix A of ref. [2] for a derivation), 4 where m W = g 2 v/2 is the W ± mass, v is the temperature-dependent Higgs expectation value, 5 and is a Debye mass [51] (for future reference we also define m 2 E1 here): For the neutral gauge field components (B 0 , A 3 0 ) the propagator is a matrix, whose form can be found in eqs. (A.22) and (A.23) of ref. [2].
With the notation introduced, the potentials appearing in eq. (4.7) read The r-independent parts, denoted somewhat formally with the argument r = 0, correspond to self-energy contributions; the r-dependent parts to exchange contributions. 6 The r-independent parts are linearly divergent, and the corresponding vacuum counterterms are defined such that lim r→∞ V i (r) = 0 at T = 0. Explicit expressions are given in appendix A. At T > 0, lim r→∞ Re[V i (r)] = 0 amounts to the Salpeter correction. As elaborated upon in sec. 4.2 and as can be deduced from eq. (4.17), in V 1 the thermal widths cancel to leading order in r ∼ 1/(M v), whereas in V 2,3,4 they represent fast reaction rates ∼ αT .

Non-degenerate situation
If λ 4 = 0 or λ 5 = 0 and v > 0, eq. (2.1) implies that different components of the inert doublet χ have different masses. A mass splitting is also induced by Standard Model radiative corrections [19]. In this situation the potentials of eqs. (4.17)-(4.19) get replaced by matrix potentials which act in the space of the field components H ± , H 0 , H0 defined in eq. (2.7). Modifying the notation slightly from eq. (2.3), we denote the mass of the neutral component H 0 by M , and the additional rest mass of the pair H i H j by ∆M ij . The kinetic masses appearing in the Schrödinger equations also depend on the pair in question, however for small but non-zero ∆M ij > ∼ α 2 M this can be considered to be a higher-order effect, and will be omitted in the following (its inclusion is trivial, by replacing the kinetic term in eq. (4.37) by a diagonal matrix containing the reduced masses). Even though eq. (3.6) contains expectation values of the type for the Schrödinger equation it is convenient to consider the opposite time ordering [46], and similarly for 1 → 2, 3, 4. The two Wightman functions are related by which is one way to see the origin of the Laplace transform in eq. (4.10). The function Π > i (ω, k) in turn agrees with the spectral function up to a trivial factor and exponentially small corrections, where n B is the Bose distribution. When the Wightman functions Π > i corresponding to the operators in eq. (2.3) are written in the basis of eq. (2.7), they have an overlap with many different "elementary" Wightman functions. The overlaps form a block-diagonal form, and can be expressed through four different "weight matrices", denoted by W i : c 2 +4(c 1 +c 3 +c 4 ) 16 Given that c 3 = c 4 (cf. eq. (2.6)), eq. (4.25) has itself a block-diagonal form.
The right-hand sides of eq. (4.7), which may be called the source terms, also turn into matrices in the basis of eq. (2.7). These matrices are diagonal, but have in some cases non-trivial coefficients, corresponding to the multiplicities of contractions: As a crosscheck, it may be noted that projecting the sources from eqs. (4.29)-(4.32) with the weights from eqs. (4.25)-(4.28) yields which indeed agrees with weighted sum over the source terms of eq. (4.7) with the normalization factors from eq. (4.9).
The potentials can be derived as explained in ref. [2], from the thermal expectation value of the time-evolution operator bracketed between states like in eqs. (4.25)-(4.28). At this point the sources are momentarily separated from each other; it is advantageous to symmetrize the state generated in this point-splitting, e.g.
Then a straightforward computation produces matrix potentials, listed in table 1. Apart from eq. (4.12), the potentials in table 1 contain the object and similarly for V ZZ and VZZ, where we have defined (4. 36) We stress that at finite temperature Z 0 does not represent a propagating mode, andZ 0 does not represent one even at zero temperature. The fields Z 0 andZ 0 simply stand for specific linear combinations originating from vertices; the diagonal modes are obtained from B 0 and A 3 0 through an orthogonal transformation parametrized by a temperature-dependent mixing angleθ, given in eq. (A.5).
The potentials of table 1 contain a real part, including the diagonal r-independent Salpeter correction, as well as an imaginary part, representing scatterings and decays of the type described by eq. (4.2). As mentioned in sec. 4.2, the inclusion of the scatterings has been demonstrated to be theoretically consistent in the case of the most attractive channel, in which case the scattering rate is a slow one. This slow rate appears in the upper diagonal block of the potential U 1 in table 1. Its role is to damp (or "decohere") oscillations between the three states appearing in this block. In the other channels, the widths represent a part of NLO corrections.
With these ingredients at hand, the thermally averaged scattering rates are obtained from matrix Schrödinger equations of the form   where the matrix F i has the same dimension as the source S i . The combination needed for eq. (3.6) becomes, in analogy with eq. (4.10), It is interesting to ask how the degenerate limit of sec. 4.3 is recovered from the equations of the current section. A simple way to do this is to recall that if a Green's function is expressed as a function of time t rather than energy E ′ , then the source terms in eqs. (4.29)-(4.32) represent initial conditions at time t = 0 [46]. To first order in interactions, we can simply act on the initial conditions with the potentials of table 1, and subsequently project the results with the weights from eqs. (4.25)-(4.28), i.e. compute 4 i=1 Tr W i U i S i . It can be verified that the terms proportional to 2c 1 δ (3) (r − r ′ ),

Limit of low temperatures
The scale hierarchy shown in eq. (4.6) breaks down as the temperature decreases: first the Debye scale gT becomes smaller than the energy scale α 2 M at which the Schrödinger description applies, and soon afterwards πT also becomes smaller than α 2 M . Moreover, assuming that mass splittings in the dark sector are ∆M ij > ∼ α 2 M , πT also becomes smaller than ∆M ij . These crossings have an important impact on the determination of i c i γ i and σ eff v at low temperatures, particularly as far as the below-threshold part (E ′ < 0) is concerned, given that the Laplace transforms in eqs. (4.10) and (4.39) exponentially enhance the contributions from the smallest energies.
It may be noted, first of all, that once the Debye scale drops below α 2 M , the dominant process responsible for the thermal interaction rate is the absorption of a thermal gauge boson (cf. fig. 2(left)) rather than scattering off Standard Model particles as is the case at higher temperatures (cf. fig. 2(right)) (cf. ref. [53] and references therein). However, this does not change the magnitude of the thermal interaction rate qualitatively. Given that the numerical effect from the low-temperature regime is modest, we have not worked out these effects quantitatively; this would pose an interesting topic for future research.
More importantly, the spectral function changes dramatically once πT < ∼ |E ′ | ∼ α 2 M . In vacuum, the spectral function vanishes for E ′ < 0 in repulsive channels, and for E ′ below the ground state energy in attractive channels. At T > 0, this is no longer the case: any "measurement" can detect non-vanishing below-threshold spectral weight, with the energy difference to the vacuum threshold supplied by a thermal fluctuation suppressed by a Boltzmann factor. This has been shown explicitly in a QCD context, both by considering the dissociation rate of bound states with pNRQCD (cf. eq. (89) of ref. [48]), and through a strict NLO computation of the process in fig. 2(left) together with the associated virtual corrections (cf. eq. (4.7) of ref. [54] after setting ω → 2M + ∆E ′ ). We have not carried out a quantitative analysis of these effects for the present system, which would again pose an interesting topic for future research, however we multiply thermal interaction rates by the Boltzmann factor θ(−E ′ )e −|E ′ |/T in order to account for the exponential suppression below threshold. This is a higher-order effect in the domain of our main interest, eq. (4.6), but imposes the correct overall magnitude to the below-threshold spectral function when πT < ∼ α 2 M .
The third effect concerns mass splittings, which are always present at least at the level ∆M ij ∼ 10 −3 M [19]. To account for them properly requires the numerical solution of the matrix equations derived in sec. 4.4. However, on the qualitative level we can profit from a corresponding solution that was worked out in sec. 7 of ref. [2]. The main finding was that as long as ∆M ij ∼ α 2 M , the shape of the spectral function does not depend noticeably on ∆M ij , however the spectral function splits into several parts, separated by the mass shifts. 7 We can work out these shifts by solving eqs. (4.37) and (4.38) at tree level but with ∆M ij = 0. 8 Denoting by ρ (0) ≡ M 3 2 θ(E ′ ) √ E ′ /(4π) the tree-level spectral function obtained with ∆M ij = 0, and using ∆M The shape stays intact because the heavier particles still contribute as virtual states and thereby generate an interaction between the lightest ones. 8 Thermal mass corrections can be omitted in this regime, given that |∆M T | ∼ α 3/2 T < ∼ α 7/2 M ≪ α 2 M .

Numerical solution and overclosure bound
Once the combination i c i γ i has been computed as a function of the temperature, either from eq. (4.10) or from eq. (4.39), the effective cross section σ eff v is obtained from eq. (3.6). Writing out the time derivative in eq. (3.1), the evolution equation reads where H is the Hubble rate. Combining this with the entropy conservation law (∂ t +3H)s = 0 as well as with the relation of time and temperature,Ṫ = −3Hs/c, where c is the heat capacity; defining a "yield parameter" through Y ≡ n/s; and denoting z ≡ M/T , we get Here m Pl is the Planck mass and e is the energy density. We insert e, c, and s from ref. [55]. Our goal is to determine a conservative overclosure bound for M . Thus, for a given M , we need a lower bound for Y . A lower bound for Y requires an upper bound for σ eff v , so that annihilations take place with maximal efficiency. As discussed in sec. 4.5, if ∆M ij ∼ α 2 M ∼ 10 −3 M , then in the non-degenerate situation the solution of the Schrödinger equation does not differ qualitatively from the degenerate limit. In fact σ eff v decreases with ∆M ij , because of the Boltzmann suppression factors ∼ e −∆M ij /T induced by the movement of the heavier particle thresholds to higher energies. Therefore, the degenerate limit sets an upper bound for σ eff v . We only depart from this approximation at very low temperatures πT < ∼ α 2 M where effects from ∆M ij start to be of order unity (cf. eqs. (5.7) and (5.8)).
For numerical evaluations, the gauge couplings g 2 1 and g 2 2 , the top Yukawa coupling h 2 t , and the scalar couplings appearing in eq. (2.1) are needed. The gauge couplings affecting the "soft" thermal physics of the static potential are evaluated at a scaleμ ≃ πT . In contrast the couplings in eqs. (2.4)-(2.6) are needed at a scaleμ ≃ 2M . We fix g 2 1 (m Z ) = 0.128, g 2 2 (m Z ) = 0.425, h 2 t (m Z ) = 0.967, λ 1 (m Z ) = 0.145, and forμ < m Z keep these unchanged. Reasonable agreement is found, in spite of the presence of Debye screening and complicated mixing patterns that appear in the thermal potentials.
Consider now σ eff v from eq. (3.6). It is convenient to express the result in a form similar to eq. (4.11), where "average Sommerfeld factors" have been defined as The Salpeter correction is given by eqs.
If the change of the threshold location were the only modification of the spectral function ρ i , 2∆M T would exactly cancel out in eq. (5.4).
As discussed in sec. 4.5, the vacuum mass differences ∆M ij become important at very low temperatures (in contrast ∆M T loses its significance there). Inserting eq. (4.40) into eq. (4.39), comparing with eq. (5.3), and setting for simplicity ∆M + = ∆M0 ≡ ∆M , the effects from ∆M can phenomenologically be included through the substitutions We adopt this recipe in the following, setting for illustration ∆M = 10 −3 M , which is parametrically in the correct range ∼ α 2 M and numerically in fair accordance with ref. [19] at λ i = 0, and also reflects the gradual increase of ∆M ≃ λ 4,5 v 2 /M with scalar self-couplings. The case ∆M = 0 is considered as an upper bound on the average Sommerfeld factors. The average Sommerfeld factors have been plotted in fig. 4. For the numerical evaluation of eq. (5.4), we have restricted the Laplace transform to the range . Given the average Sommerfeld factors, we can insert eq. (5.3) into eq. (5.2) and integrate the latter equation for Y (z). Examples of solutions are shown in fig. 5. We have compared with the linearized version of this equation (cf. eq. (3.4)), obtained by setting Y 2 − Y 2 eq → 2Y eq (Y − Y eq ). It is observed how the initial departure from equilibrium is well described by both forms, however afterwards the Lee-Weinberg from of eq. (5.2) leads to a substantial depletion of the dark matter abundance.
As can be deduced from fig. 5, Y eq has become exponentially small by the time that z ∼ 40. In the absence of Y eq , eq. (5.2) can be integrated into The regime z > ∼ 40 can easily reduce the dark matter abundance by a factor 2...3. We choose z final = 10 4 so that the contribution from late times is typically at the percent level. Note that weak interactions are faster than the Hubble rate down to T ≃ 10 MeV, so we may assume the dark matter particles to be kinetically equilibrated in the whole z range. It should however be noted that, taken literally, the growing Sommerfeld factorS 1,eff in fig. 4 compromises the convergence of eq. (5.9) at large z. At the same time, at low temperatures kinetic and chemical equilibrium is gradually lost in the dark sector, and the bound-state thermal abundance is presumably no longer available as an efficient annihilation channel once πT ≪ α 2 M . The value z final = 10 4 represents a phenomenological compromise where the numerical effect from large z is small, yet the physics assumptions that went into the thermal analysis should still be intact. It would be interesting to understand the physics of this regime more precisely (cf. also the comments in secs. 4.5 and 6).
Eventually the heavier dark matter particles decay into the lightest one, so that the final yield is Y phys = Y (z final ). The energy density carried by the lightest ones today is ρ dm (T 0 ) = M Y phys s(T 0 ), and the energy fraction is Ω dm (T 0 ) = M Y phys s(T 0 )/ρ cr (T 0 ), where ρ cr is the current critical energy density. Inserting from ref. [58] s(T 0 ) = 2 891/cm 3  which can be compared with the observed value Ω dm h 2 obs = 0.1186(20) [59]. Results are plotted in fig. 6; a discussion is deferred to the first paragraph of sec. 6.

Conclusions and outlook
The purpose of this paper has been to illustrate and refine the general formalism of ref. [2], by applying it to a simple yet phenomenologically viable dark matter computation. After the inclusion of thermal effects, such as the Salpeter correction to dark matter masses, the modification of the Sommerfeld effect through Debye screening, and thermal interaction rates, we find a conservative upper bound for the mass of the lightest dark matter particle within the inert doublet model (IDM), as a function of quartic scalar couplings. As a reference, we note that for vanishing quartic couplings values M < ∼ 535 ± 9 GeV can typically be found in literature (cf. e.g. refs. [25,36]), and that for this case we get M < ∼ 519 ± 4 GeV by using free spectral functions (cf. fig. 6(right)). Switching on the thermally modified Sommerfeld factors, Salpeter corrections, and thermal interaction rates, the bound increases to M < ∼ 523± 5 GeV for ∆M = 10 −3 M , and to M < ∼ 562 ± 5 GeV for the extreme case ∆M/M → 0. For the maximal quartic couplings considered, λ 3 (2M ) = λ 4 (2M ) = λ 5 (2M ) = π, we obtain M < ∼ 10.6 ± 0.1 TeV with free spectral functions; M < ∼ 11.1 ± 0.1 TeV for ∆M = 10 −3 M ; and M < ∼ 12.5 ± 0.1 TeV for ∆M/M → 0. The uncertainties cited here originate from the observed value of the dark matter relic density [59].
In the high-mass regime the system displays a non-trivial bound-state spectrum at low temperatures (cf. fig. 6(left)), which leads to large Sommerfeld factors at large z (cf. fig. 4). This results in efficient annihilation, and helps to push up the upper bound for M . We stress that the bound-state spectrum is easily addressed within our formalism, since the known Hard Thermal Loop resummed thermal interaction rate (reflecting the processes in fig. 2(right)) eliminates the need for complicated bound-state production and dissociation rate computations. At very low temperatures, T ≪ α 3/2 M , other processes contribute as well (cf. fig. 2(left)), however these have also been studied in the QCD context (cf. refs. [53,54] and references therein), and the same techniques could conceivably be generalized to cosmology. Once πT ≪ α 2 M , there is gradual departure from kinetic and chemical equilibrium in the dark sector, whose study represents a complicated but interesting open problem. 10 Once the collider lower bound exceeds the cosmological upper bound of fig. 6(right), IDM is firmly excluded as a model, independently of astrophysical uncertainties related to the local dark matter distribution. In practice, accepting modest astrophysical assumptions, direct and indirect non-detection constraints permit to set more stringent bounds than the overclosure one (cf. e.g. refs. [35,36] and references therein).
One weakness of the IDM is that the quartic scalar couplings can be varied in a broad range, which has a significant effect on the overclosure bound (cf. fig. 6(right)). The quartic couplings also influence mass splittings, resulting in a non-trivial multidimensional parameter dependence. If the couplings are large, their effects should be resummed. For instance the scalar couplings affect the thermal corrections to dark matter masses; in contrast to the Salpeter correction in eq. (5.5), these effects are power-suppressed, ∆M T ≃ (2λ 3 + λ 4 )T 2 /(24M ). In addition, at T < ∼ 160 GeV, the Higgs mechanism (v > 0) generates cubic scalar couplings which lead to additional terms in the static potentials (cf. e.g. ref. [36]). In the present investigation we resummed only effects from gauge couplings, which are not suppressed by T /M or v/M and are therefore expected to generically give the dominant contributions.

Acknowledgements
This work was supported by the Swiss National Science Foundation (SNF) under grant 200020-168988. S.B. thanks Germano Nardini and Lewis Tunstall for helpful discussions.
(A. 15) Defining α 1 ≡ (3g 2 2 + g 2 1 )/(16π), α 2 ≡ (g 2 2 − g 2 1 )/(16π) and α 3,4 ≡ (g 2 2 + g 2 1 )/(16π), the corresponding Sommerfeld factors read [10] S 1 = X 1 1 − e −X 1 , S 2,3,4 = X 2,3,4 e X 2,3,4 − 1 , (A. 16) where X i ≡ πα i /v and v parametrizes E ′ from eq. (4.10) as E ′ = 2∆M T + M v 2 . We note that eqs. (A.8)-(A.13) are based on evaluating gauge field self-energies in the Hard Thermal Loop approximation. This is justified as long as the particles with which gauge fields interact are ultrarelativistic, i.e. with masses m ≪ πT . If m > ∼ πT , the self-energies take a more complicated form (cf. appendix A of ref. [2] for the full 1-loop self-energy matrix of the neutral components A 3 0 , B 0 ), and thermal modifications cannot be captured by the two Debye mass parameters m 2 E1 and m 2 E2 . Nevertheless, it is possible to identify the light-fermion contribution to the Debye masses. If we consider vanishing spatial momentum; model top and bottom quarks by a common "fermionic" mass m f ; and model W ± , Z 0 and Higgs bosons by a common "gauge" mass m g ; then eq. (A.6) of ref. [52] shows that terms mixing A 3 0 and B 0 drop out, and we may replace eq. (4.16) with Here the fermionic and bosonic susceptibilities read where n F and n B are the Fermi and Bose distributions, respectively. We have adopted eq. (A.17) for modelling the low-temperature regime, inserting m f ≃ (m t m b ) 1/2 and m g ≃ (m Z m 2 W m φ ) 1/4 , but stress that this represents a purely phenomenological recipe within the complicated temperature interval m b < ∼ πT < ∼ m t .