Re-derived overclosure bound for the inert doublet model

TL;DR

This paper refines the overclosure bound for the inert doublet dark matter model by incorporating thermal effects and higher-order corrections, revealing a modest increase in the bound across a broad temperature range.

Contribution

It introduces a formalism that accounts for thermal effects and resummed corrections in calculating dark matter overclosure bounds for the inert doublet model.

Findings

Resummed corrections increase the overclosure bound by 1-18%.

Thermal effects are significant in late annihilation stages.

The analysis covers a broad temperature range, emphasizing late-stage effects.

Abstract

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.

Tables

Table 1. Table 1: The “potentials” 𝒰 i subscript 𝒰 𝑖 \mathcal{U}_{i} appearing in eq. ( 4.72 ). In general the potentials contain both a real part, as well as an imaginary part representing thermal scatterings (cf. eq. ( 4.15 ) and sec.  4.2 ).

${𝒰}_1\left(r\right)$
$⟨H_{+}{H}_{-}$	$⟨H_0{H}_0$	$⟨H_{\overline{0}}{H}_{\overline{0}}$	$⟨i{H}_0{H}_{\overline{0}}$
${𝒱}_{\overline{Z}\overline{Z}}\left(0\right)+2{𝒱}_{\mathrm{WW}}\left(0\right)-{𝒱}_{\overline{Z}\overline{Z}}\left(r\right)$	$-{𝒱}_{\mathrm{WW}}\left(r\right)$	$-{𝒱}_{\mathrm{WW}}\left(r\right)$	0	$H_{+}^{†}{H}_{-}^{†}⟩$
$-2{𝒱}_{\mathrm{WW}}\left(r\right)$	${𝒱}_{\mathrm{ZZ}}\left(0\right)+2{𝒱}_{\mathrm{WW}}\left(0\right)$	$-{𝒱}_{\mathrm{ZZ}}\left(r\right)$	0	$H_0^{†}{H}_0^{†}⟩$
$-2{𝒱}_{\mathrm{WW}}\left(r\right)$	$-{𝒱}_{\mathrm{ZZ}}\left(r\right)$	${𝒱}_{\mathrm{ZZ}}\left(0\right)+2{𝒱}_{\mathrm{WW}}\left(0\right)$	0	$H_{\overline{0}}^{†}{H}_{\overline{0}}^{†}⟩$
0	0	0	${𝒱}_{\mathrm{ZZ}}\left(0\right)+2{𝒱}_{\mathrm{WW}}\left(0\right)+{𝒱}_{\mathrm{ZZ}}\left(r\right)$	$-i{H}_0^{†}{H}_{\overline{0}}^{†}⟩$

Equations

$${\mathcal{L}}_{\chi }$$

$$\chi=\frac{1}{\sqrt{2M}}\Bigl{(}Ce^{-iMt}+D^{\dagger}e^{iMt}\Bigr{)}\;,\quad\chi^{\dagger}=\frac{1}{\sqrt{2M}}\Bigl{(}De^{-iMt}+C^{\dagger}e^{iMt}\Bigr{)}\;,$$

$$\delta\mathcal{L}_{\mbox{\scriptsize abs}}=i\,\bigl{(}c_{1}\,\underbrace{C^{\dagger}_{p}D^{\dagger}_{p}D_{q}C_{q}}_{\equiv\;O_{1}}\,+\,c_{2}\,\underbrace{C^{\dagger}_{p}T^{a}_{pq}D^{\dagger}_{q}\,D_{r}T^{a}_{rs}C_{s}}_{\equiv\;O_{2}}\,+\,c_{3}\,\underbrace{D^{\dagger}_{p}D^{\dagger}_{q}D_{p}D_{q}}_{\equiv\;O_{3}}\,+\,c_{4}\,\underbrace{C^{\dagger}_{p}C^{\dagger}_{q}C_{p}C_{q}}_{\equiv\;O_{4}}\;\bigr{)}\;.$$

$$c_1$$

$$c_2$$

$$c_3$$

$$C=\left(\begin{array}[]{c}H_{+}\\[2.84526pt] {\displaystyle\frac{H_{0}-iH_{\bar{0}}}{\sqrt{2}}}\end{array}\right)\;,\quad D=\left(\begin{array}[]{c}H_{-}\\[2.84526pt] {\displaystyle\frac{H_{0}+iH_{\bar{0}}}{\sqrt{2}}}\end{array}\right)\;.$$

$$\dot{n}\;=\;-\langle\sigma_{\mbox{\scriptsize eff}}\,v\rangle\,\bigl{(}n^{2}-n_{\mbox{\scriptsize eq}}^{2}\bigr{)}\;,$$

$$⟨{\sigma }_{\text{mbox}eff}v⟩=\sum _{i,j}\frac{⟨{\sigma }_{ij}{v}_{ij}⟩n_i^{\text{mbox}eq}{n}_j^{\text{mbox}eq}}{n_{\text{mbox}eq}^2}$$

$$n_{\mbox{\scriptsize eq}}\;\approx\;\int_{{\bf k}}\Bigl{(}e^{-E_{+}/T}+e^{-E_{-}/T}+e^{-E_{0}/T}+e^{-E_{\bar{0}}/T}\Bigr{)}\;=\;\sum_{i=\pm,0,\bar{0}}\frac{TM_{H_{i}}^{2}}{2\pi^{2}}K_{2}\biggl{(}\frac{M_{H_{i}}}{T}\biggr{)}\;,$$

$$\dot{n}\;=\;-\Gamma_{\mbox{\scriptsize chem}}\bigl{(}n-n_{\mbox{\scriptsize eq}}\bigr{)}+{\mathcal{O}}\bigl{(}n-n_{\mbox{\scriptsize eq}}\bigr{)}^{2}\;,$$

$$⟨{\sigma }_{\text{mbox}eff}v⟩\equiv \frac{{\Gamma }_{\text{mbox}chem}}{2n_{\text{mbox}eq}}.$$

$$\langle\sigma_{\mbox{\scriptsize eff}}\,v\rangle=\frac{4}{n_{\mbox{\scriptsize eq}}^{2}}\sum_{i=1}^{4}c_{i}\gamma_{i}\;,\quad\gamma_{i}\;\equiv\;\bigl{\langle}O_{i}\bigr{\rangle}\;.$$

$${\chi }_i{\chi }_j\leftrightarrow X{X}^{\prime },$$

$${\chi }_{i}X\leftrightarrow {\chi }_j{X}^{\prime },{\chi }_i\leftrightarrow {\chi }_{j}X{X}^{\prime }.$$

$$({\chi }_i{\chi }_k{)}_{\text{mbox}open}X\leftrightarrow ({\chi }_j{\chi }_l{)}_{\text{mbox}bound}{X}^{\prime },({\chi }_i{\chi }_k{)}_{\text{mbox}open}\leftrightarrow ({\chi }_j{\chi }_l{)}_{\text{mbox}bound}X{X}^{\prime },$$

$$({\chi }_i{\chi }_k{)}_{\text{mbox}open}X\leftrightarrow ({\chi }_j{\chi }_l{)}_{\text{mbox}open}{X}^{\prime },({\chi }_i{\chi }_k{)}_{\text{mbox}open}\leftrightarrow ({\chi }_j{\chi }_l{)}_{\text{mbox}open}X{X}^{\prime },$$

$${\chi }_i{\chi }_j(\text{mbox}virtualX)\leftrightarrow X^{\prime }{X}^{\prime \prime },$$

$${\alpha }^{2}M\ll gT,\alpha M,\pi T\ll M,$$

$$\displaystyle\biggl{[}-\frac{\nabla_{r}^{2}}{M}+\mathcal{V}_{i}(r)-E^{\prime}\biggr{]}G_{i}(E^{\prime};{\bf r},{\bf r^{\prime}})$$

$$\underset{\mathbf{r},{\mathbf{r}}^{\prime }\to 0}{\lim }\text{mbox}Im{G}_i(E^{\prime };\mathbf{r},{\mathbf{r}}^{\prime })$$

$$N_1=2,N_2=\frac{3}{2},N_3=N_4=6.$$

$${\gamma }_i$$

$$⟨{\sigma }_{\text{mbox}eff}v{⟩}^{(0)}=\frac{c_1}{2}+\frac{3c_2}{8}+\frac{3(c_3+c_4)}{2}.$$

$$\displaystyle\mathcal{V}_{\mbox{\tiny\rm{${W}{W}$}}}(r)$$

$$\displaystyle\mathcal{V}_{\mbox{\tiny\rm{${A}{A}$}}}(r)$$

$$\displaystyle\mathcal{V}_{\mbox{\tiny\rm{${B}{B}$}}}(r)$$

$$\displaystyle i\bigl{\langle}W^{+}_{0}W^{-}_{0}\bigr{\rangle}_{\mbox{\tiny\rm{T}}}(0,k)$$

$$m^{2}_{\mbox{\tiny\rm{E1}}}\;\equiv\;\Bigl{(}{\frac{n_{\mbox{\tiny\rm{$S$}}}}{6}}+\frac{5n_{\mbox{\tiny\rm{$G$}}}}{9}\Bigr{)}g_{1}^{2}T^{2}\;,\quad m^{2}_{\mbox{\tiny\rm{E2}}}\;\equiv\;\Bigl{(}{\frac{2}{3}}+{\frac{n_{\mbox{\tiny\rm{$S$}}}}{6}}+\frac{n_{\mbox{\tiny\rm{$G$}}}}{3}\Bigr{)}g_{2}^{2}T^{2}\;,\quad n_{\mbox{\tiny\rm{$S$}}}\equiv 1\;,\quad n_{\mbox{\tiny\rm{$G$}}}\equiv 3\;.$$

$${\mathcal{V}}_1(r)$$

$${\mathcal{V}}_2(r)$$

$${\mathcal{V}}_{3,4}(r)$$

$${\gamma }_i$$

$${\Pi }_1^{<}(\omega ,\mathbf{k})$$

Full text

July 2017

{centering}

**Re-derived overclosure bound for the inert doublet model **

S. Biondini and M. Laine

*AEC, Institute for Theoretical Physics, University of Bern,

Sidlerstrasse 5, CH-3012 Bern, Switzerland

Abstract

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\sim 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.

1 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\mathop{\raise 1.29167pt\hbox{>\kern-7.5pt\raise-4.73611pt\hbox{\sim}}}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],111Other discussions of thermal effects relevant for heavy particles can be found e.g. in refs. [3, 4, 5, 6, 7, 8].

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\sim 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^{\pm},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 co-annihilation 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 $\sim\nabla^{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. [19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41] 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).

2 4-particle operators

In the inert doublet model (IDM), the Standard Model is supplemented by an additional Higgs doublet, $\chi$, which does not couple to fermions because of an unbroken discrete Z(2) symmetry. Denoting by $\phi$ the Standard Model Higgs doublet and by $D_{\mu}$ the corresponding covariant derivative, the Standard Model Lagrangian is modified by the additional terms

[TABLE]

The notation $\lambda_{1}$ is reserved for the Standard Model Higgs self-coupling, $\delta\mathcal{L}_{\mbox{\tiny\rm{SM}}}=-\lambda_{1}(\phi^{\dagger}\phi)^{2}$.

If the mass scale $M$ is much larger than the electroweak scale, M\gg m_{\mbox{\tiny\rm{W}}}, the $\chi$ particles annihilate efficiently into the Standard Model ones. The annihilations that we are interested in happen in the temperature range $T\sim M/20...M/10^{4}$, in which the average velocity is $v\sim\sqrt{T/M}\ll 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

[TABLE]

then at leading order in $1/M^{2}$ there are four “absorptive” operators that play a role:222As 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].

[TABLE]

Here sums over the isospin components $p,q,r,s\in\{1,2\}$ are implied, and $T^{a}\equiv\sigma^{a}/2$, where $\sigma^{a}$ are the Pauli matrices.

We have computed the coefficients $c_{1},...,c_{4}$ in eq. (2.3) in general $R_{\xi}$ gauges at leading non-trivial order, verifying their gauge independence for \xi<M^{2}/m_{\mbox{\tiny\rm{Z}}}^{2}:333For \xi\gg M^{2}/m_{\mbox{\tiny\rm{Z}}}^{2}\gg 1 the results change qualitatively and therefore unitary gauge is not viable.

[TABLE]

Here $g_{1}$ and $g_{2}$ are the U${}_{\mbox{\tiny\rm{Y}}}$(1) and SU${}_{\mbox{\tiny\rm{L}}}$(2) gauge couplings, respectively. The couplings should be evaluated at a renormalization scale $\sim 2M$. The same values of the coefficients can be extracted from ref. [25].

If $\lambda_{4}\neq 0$ or $\lambda_{5}\neq 0$, or if Standard Model radiative corrections are considered, different components of $\chi$ are non-degenerate in mass. In this case the doublets $C$ and $D$ can be written as

[TABLE]

The operators in eq. (2.3) split into a $10\times 10$ matrix in the field space of eq. (2.7), which is given (with a slightly different notation) in eqs. (4.31)–(4.46) below.

3 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\equiv\sum_{i=\pm,0,\bar{0}}n_{i}$. Within a Boltzmann approach, ref. [15] established that $n$ evolves according to the Lee-Weinberg equation [43, 44],

[TABLE]

where $\dot{n}$ is the covariant time derivative in an expanding background, and

[TABLE]

is an effective cross section for $2\to 2$ annihilations from the dark sector. In our case the total equilibrium number density reads, at tree-level,

[TABLE]

where $\int_{{\bf k}}\equiv\int\frac{{\rm d}^{3}{\bf k}}{(2\pi)^{3}}$, $E_{i}\equiv\sqrt{k^{2}+M_{H_{i}}^{2}}$ with $k\equiv|{\bf k}|$, 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 \Delta M_{\mbox{\tiny\rm{T}}}\sim-\alpha^{3/2}T<0, where $\alpha$ 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]

[TABLE]

where $\Gamma_{\mbox{\scriptsize 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

[TABLE]

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 $\Gamma_{\mbox{\scriptsize 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), $\Gamma_{\mbox{\scriptsize chem}}$ is to first order proportional to the thermal expectation value of $\delta\mathcal{L}_{\mbox{\scriptsize abs}}$ [9]. Inserting the proportionality coefficient and expressing the result through eq. (3.5), we obtain

[TABLE]

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 $\gamma_{i}$ is exponentially suppressed by $\sim e^{-2M/T}$, with the Boltzmann factor cancelling against that from $n_{\mbox{\scriptsize eq}}^{2}$.

Eq. (3.6) represents a generalization of eq. (3.2). The matrix structure of $\sigma_{ij}$ corresponds to matrix-like Schrödinger equations satisfied by the wave functions of the annihilating pair (cf. table 1), and the weights $n^{\mbox{\scriptsize eq}}_{i}$ in eq. (3.2) correspond to threshold locations in the Laplace transform in eq. (4.74). 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.

4 Schrödinger description

4.1 General goal and physical interpretation

In the notation of ref. [15], the cross sections in eq. (3.2) describe the processes

[TABLE]

where $X,X^{\prime}$ 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 $\Gamma\sim\frac{\alpha^{2}}{M^{2}}\int_{{\bf k}}e^{-E/T}$. However, the $\chi_{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]:

[TABLE]

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_{\mbox{\scriptsize eq}}^{i}$ into eq. (3.2). If there are bound states in the dark sector, further “fast” processes can be added, notably

[TABLE]

where we assume that the binding energy is small, \Delta E\sim\alpha^{2}M\mathop{\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}}\pi T. Of course the same reactions are also present without bound states,

[TABLE]

and can change the annihilating pair into a different gauge or spin state. In addition, processes with virtual $X$ exchange are important,

[TABLE]

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.

4.2 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 $\Delta E\sim\alpha^{2}M$ provided that (cf. e.g. refs. [46, 47, 48, 49])

[TABLE]

where $\alpha\sim g^{2}/(4\pi)$. 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 $\sim{\bf r}\cdot g{\bf E}$. 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 ${\bf 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, $\Gamma\sim\alpha^{2}T^{3}r^{2}$. Therefore, $\Gamma$ 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 $\chi$ pairs in gauge non-singlet channels, are of order $\alpha T$. This is a fast rate, rapidly decohering the phase of the wave function and justifying a classical Boltzmann description. At the same time, it is not clear whether such a rate can be consistently included in a Schrödinger equation: if T\mathop{\raise 1.29167pt\hbox{>\kern-7.5pt\raise-4.73611pt\hbox{\sim}}}\alpha M, $\Gamma\sim\alpha T$ modifies the spectral function in the energy range $\Delta E\sim\alpha^{2}M$ by an effect of ${\mathcal{O}}(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.

4.3 Degenerate limit

We start by considering the degenerate limit, i.e. $M\equiv M_{0}=M_{\bar{0}}=M_{\pm}$. Each of the expectation values in eq. (3.6) can be expressed as a Laplace transform of a spectral function, denoted by $\rho_{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]:

[TABLE]

where $N_{i}$ is a normalization factor giving the number of contractions related to $O_{i}$:

[TABLE]

In center-of-mass coordinates the Laplace transform reads

[TABLE]

where $M\gg\Lambda\gg\alpha^{2}M$ is a cutoff restricting the average to the non-relativistic regime. According to eq. (3.6), the physical result is $\sum_{i=1}^{4}c_{i}\gamma_{i}$, with $c_{i}$ given in eqs. (2.4)–(2.6).

In the free limit, $\mathcal{V}_{i}\to 0$, the spectral function from eqs. (4.7) and (4.8) reads $\rho^{(0)}_{i}(E^{\prime})=N_{i}M^{{\frac{3}{2}}}\theta(E^{\prime})\sqrt{E^{\prime}}/(4\pi)$. Carrying out the Laplace transform in eq. (4.10), inserting n^{(0)}_{\mbox{\scriptsize eq}}=4\bigl{(}\frac{MT}{2\pi}\bigr{)}^{{\frac{3}{2}}}e^{-M/T} from eq. (3.3), and plugging into eq. (3.6), we obtain the value of $\langle\sigma_{\mbox{\scriptsize eff}}\,v\rangle$ for a degenerate system and to leading order in $1/M^{2}$ and $\alpha$:

[TABLE]

In order to go beyond eq. (4.11), we include the potentials $\mathcal{V}_{i}$ for the various channels in eq. (4.7). It is helpful to introduce the notation

[TABLE]

where $\langle...\rangle_{\mbox{\tiny\rm{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),444For hard momenta k\gg m_{\mbox{\tiny\rm{\widetilde{W}}}} only the massless part $1/k^{2}$ is important. The full form is needed for correctly estimating the contribution of soft near-threshold momenta to the annihilation cross section. The soft momenta become increasingly important as the temperature decreases.

[TABLE]

where m_{\mbox{\tiny\rm{W}}}=g_{2}v/2 is the $W^{\pm}$ mass, $v$ is the temperature-dependent Higgs expectation value,555Even though carrying the same symbol, $v$ should not be confused with the non-relativistic velocity appearing e.g. in eq. (3.1).

and m_{\mbox{\tiny\rm{\widetilde{W}}}}^{2}\equiv m_{\mbox{\tiny\rm{W}}}^{2}+m_{\mbox{\tiny\rm{E2}}}^{2}, where $m_{\mbox{\tiny\rm{E2}}}^{2}$ is a Debye mass [51] (for future reference we also define $m_{\mbox{\tiny\rm{E1}}}^{2}$ here):

[TABLE]

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

[TABLE]

The $r$-independent parts, denoted somewhat formally with the argument $r=0$, correspond to self-energy contributions; the $r$-dependent parts to exchange contributions.666The $r$-dependent parts vanish at $r\to\infty$, so that \lim_{r\to\infty}\mathcal{V}_{i}(r)=2\mathcal{V}_{\mbox{\tiny\rm{{W}{W}}}}(0)+\mathcal{V}_{\mbox{\tiny\rm{{A}{A}}}}(0)+\mathcal{V}_{\mbox{\tiny\rm{{B}{B}}}}(0).

The $r$-independent parts are linearly divergent, and the corresponding vacuum counterterms are defined such that $\lim_{r\to\infty}\mathcal{V}_{i}(r)=0$ at $T=0$. Explicit expressions are given in appendix A. At $T>0$, $\lim_{r\to\infty}\mathop{\mbox{Re}}[\mathcal{V}_{i}(r)]\neq 0$ amounts to the Salpeter correction. As elaborated upon in sec. 4.2 and as can be deduced from eq. (4.17), in $\mathcal{V}_{1}$ the thermal widths cancel to leading order in $r\sim 1/(Mv)$, whereas in $\mathcal{V}_{2,3,4}$ they represent fast reaction rates $\sim\alpha T$.

4.4 Non-degenerate situation

If $\lambda_{4}\neq 0$ or $\lambda_{5}\neq 0$ and $v>0$, eq. (2.1) implies that different components of the inert doublet $\chi$ 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_{\pm}$, $H_{0}$, $H_{\bar{0}}$ 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 $\Delta M_{ij}$. The kinetic masses appearing in the Schrödinger equations also depend on the pair in question, however for small but non-zero \Delta M_{ij}\mathop{\raise 1.29167pt\hbox{>\kern-7.5pt\raise-4.73611pt\hbox{\sim}}}\alpha^{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.72) by a diagonal matrix containing the reduced masses).

Even though eq. (3.6) contains expectation values of the type

[TABLE]

for the Schrödinger equation it is convenient to consider the opposite time ordering [46],

[TABLE]

and similarly for $1\to 2,3,4$. The two Wightman functions are related by

[TABLE]

which is one way to see the origin of the Laplace transform in eq. (4.10). The function $\Pi^{>}_{i}(\omega,{\bf k})$ in turn agrees with the spectral function up to a trivial factor and exponentially small corrections,

[TABLE]

where $n_{\mbox{\tiny\rm{B{}}}}$ is the Bose distribution.

When the Wightman functions $\Pi^{>}_{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 $\mathcal{W}_{i}$:

[TABLE]

Given that $c_{3}=c_{4}$ (cf. eq. (2.6)), eq. (4.31) 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:

[TABLE]

As a crosscheck, it may be noted that projecting the sources from eqs. (4.53)–(4.67) with the weights from eqs. (4.31)–(4.46) yields

[TABLE]

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.31)–(4.46). 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.

[TABLE]

Then a straightforward computation produces matrix potentials, listed in table 1.

Apart from eq. (4.12), the potentials in table 1 contain the object

[TABLE]

and similarly for \mathcal{V}_{\mbox{\tiny\rm{{Z}{Z}}}} and \mathcal{V}_{\mbox{\tiny\rm{\bar{Z}\bar{Z}}}}, where we have defined

[TABLE]

We stress that at finite temperature $Z_{0}$ does not represent a propagating mode, and $\bar{Z}_{0}$ does not represent one even at zero temperature. The fields $Z_{0}$ and $\bar{Z}_{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 $\tilde{\theta}$, 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 $\mathcal{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

[TABLE]

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),

[TABLE]

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^{\prime}$, then the source terms in eqs. (4.53)–(4.67) 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.31)–(4.46), i.e. compute \sum_{i=1}^{4}{\rm Tr\,}\bigl{[}\mathcal{W}_{i}\,\mathcal{U}_{i}S_{i}\bigr{]}. It can be verified that the terms proportional to $2c_{1}\,\delta^{(3)}({\bf r-r^{\prime}})$, $\frac{3c_{2}}{2}\,\delta^{(3)}({\bf r-r^{\prime}})$, and $6\,(c_{3}+c_{4})\,\delta^{(3)}({\bf r-r^{\prime}})$ reproduce the potentials from eqs. (4.17), (4.18) and (4.19), respectively.

4.5 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 $\alpha^{2}M$ at which the Schrödinger description applies, and soon afterwards $\pi T$ also becomes smaller than $\alpha^{2}M$. Moreover, assuming that mass splittings in the dark sector are \Delta M_{ij}\mathop{\raise 1.29167pt\hbox{>\kern-7.5pt\raise-4.73611pt\hbox{\sim}}}\alpha^{2}M, $\pi T$ also becomes smaller than $\Delta M_{ij}$. These crossings have an important impact on the determination of $\sum_{i}c_{i}\gamma_{i}$ and $\langle\sigma_{\mbox{\scriptsize eff}}\,v\rangle$ at low temperatures, particularly as far as the below-threshold part ($E^{\prime}<0$) is concerned, given that the Laplace transforms in eqs. (4.10) and (4.74) exponentially enhance the contributions from the smallest energies.

It may be noted, first of all, that once the Debye scale drops below $\alpha^{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 \pi T\mathop{\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}}|E^{\prime}|\sim\alpha^{2}M. In vacuum, the spectral function vanishes for $E^{\prime}<0$ in repulsive channels, and for $E^{\prime}$ 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 $\omega\to 2M+\Delta E^{\prime}$). 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 $\theta(-E^{\prime})e^{-|E^{\prime}|/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 \pi T\mathop{\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}}\alpha^{2}M.

The third effect concerns mass splittings, which are always present at least at the level $\Delta M_{ij}\sim 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 $\Delta M_{ij}\sim\alpha^{2}M$, the shape of the spectral function does not depend noticeably on $\Delta M_{ij}$, however the spectral function splits into several parts, separated by the mass shifts.777The shape stays intact because the heavier particles still contribute as virtual states and thereby generate an interaction between the lightest ones.

We can work out these shifts by solving eqs. (4.72) and (4.73) at tree level but with $\Delta M_{ij}\neq 0$.888Thermal mass corrections can be omitted in this regime, given that |\Delta M_{\mbox{\tiny\rm{T}}}|\sim\alpha^{3/2}T\mathop{\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}}\alpha^{7/2}M\ll\alpha^{2}M.

Denoting by $\rho^{(0)}\equiv M^{{\frac{3}{2}}}\theta(E^{\prime})\sqrt{E^{\prime}}/(4\pi)$ the tree-level spectral function obtained with $\Delta M_{ij}=0$, and using $\Delta M_{+}=\Delta M_{-}$ (here $\Delta M_{i_{1}..i_{n}}\equiv M_{i_{1}}+...+M_{i_{n}}-nM$), we find

[TABLE]

Inserting this into eq. (4.74), the contributions of the shifted thresholds get suppressed by $e^{-\Delta M_{ij}/T}$ just like in eq. (3.2); a practical implementation is shown in eqs. (5.7) and (5.8).

5 Numerical solution and overclosure bound

Once the combination $\sum_{i}c_{i}\gamma_{i}$ has been computed as a function of the temperature, either from eq. (4.10) or from eq. (4.74), the effective cross section $\langle\sigma_{\mbox{\scriptsize eff}}v\rangle$ is obtained from eq. (3.6). Writing out the time derivative in eq. (3.1), the evolution equation reads

[TABLE]

where $H$ is the Hubble rate. Combining this with the entropy conservation law $(\partial_{t}+3H)s=0$ as well as with the relation of time and temperature, $\dot{T}=-{3Hs}/{c}$, where $c$ is the heat capacity; defining a “yield parameter” through $Y\equiv n/s$; and denoting $z\equiv M/T$, we get

[TABLE]

Here $m_{\mbox{\scriptsize 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 $\langle\sigma_{\mbox{\scriptsize eff}}\,v\rangle$, so that annihilations take place with maximal efficiency. As discussed in sec. 4.5, if $\Delta M_{ij}\sim\alpha^{2}M\sim 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 $\langle\sigma_{\mbox{\scriptsize eff}}\,v\rangle$ decreases with $\Delta M_{ij}$, because of the Boltzmann suppression factors $\sim e^{-\Delta M_{ij}/T}$ induced by the movement of the heavier particle thresholds to higher energies. Therefore, the degenerate limit sets an upper bound for $\langle\sigma_{\mbox{\scriptsize eff}}\,v\rangle$. We only depart from this approximation at very low temperatures \pi T\mathop{\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}}\alpha^{2}M where effects from $\Delta M_{ij}$ start to be of order unity (cf. eqs. (5.7) and (5.8)).

For numerical evaluations, the gauge couplings $g_{1}^{2}$ and $g_{2}^{2}$, the top Yukawa coupling $h_{t}^{2}$, 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 $\bar{\mu}\simeq\pi T$. In contrast the couplings in eqs. (2.4)–(2.6) are needed at a scale $\bar{\mu}\simeq 2M$. We fix g_{1}^{2}(m_{\mbox{\tiny\rm{Z}}})=0.128, g_{2}^{2}(m_{\mbox{\tiny\rm{Z}}})=0.425, h_{t}^{2}(m_{\mbox{\tiny\rm{Z}}})=0.967, \lambda_{1}(m_{\mbox{\tiny\rm{Z}}})=0.145, and for \bar{\mu}<m_{\mbox{\tiny\rm{Z}}} keep these unchanged. For m_{\mbox{\tiny\rm{Z}}}<\bar{\mu}<M, the couplings are evolved like in the Standard Model, e.g. $g_{1}^{2}(\bar{\mu})\approx 48\pi^{2}/[41\ln(\Lambda_{1}/\bar{\mu})]$ and $g_{2}^{2}(\bar{\mu})\approx 48\pi^{2}/[19\ln(\bar{\mu}/\Lambda_{2})]$. For $\bar{\mu}>M$ we switch to the IDM evolution [56], $g_{1}^{2}(\bar{\mu})\approx 48\pi^{2}/[42\ln(\Lambda^{\prime}_{1}/\bar{\mu})]$ and $g_{2}^{2}(\bar{\mu})\approx 48\pi^{2}/[18\ln(\bar{\mu}/\Lambda^{\prime}_{2})]$.

Examples of spectral functions from eq. (4.8), for the three potentials from eqs. (4.17)–(4.19), are shown in fig. 3 for $M=4$ TeV, $T=M/20$.999For the numerical solution we employ the same method as in ref. [2], originally introduced in ref. [57].

The results are compared with massless Sommerfeld factors from eq. (A.16), shifted by the Salpeter correction in eq. (5.5). Reasonable agreement is found, in spite of the presence of Debye screening and complicated mixing patterns that appear in the thermal potentials.

Consider now $\langle\sigma_{\mbox{\scriptsize eff}}\,v\rangle$ from eq. (3.6). It is convenient to express the result in a form similar to eq. (4.11),

[TABLE]

where “average Sommerfeld factors” have been defined as

[TABLE]

The Salpeter correction is given by eqs. (4.17)–(4.19), (A.9), (A.11), and (A.13),

[TABLE]

Its appearance in eq. (5.4) originates from the fact that $1/n_{\mbox{\scriptsize eq}}^{2}$ in eq. (3.6) gets changed,

[TABLE]

If the change of the threshold location were the only modification of the spectral function $\rho_{i}$, 2\Delta M_{\mbox{\tiny\rm{T}}} would exactly cancel out in eq. (5.4).

As discussed in sec. 4.5, the vacuum mass differences $\Delta M_{ij}$ become important at very low temperatures (in contrast \Delta M_{\mbox{\tiny\rm{T}}} loses its significance there). Inserting eq. (LABEL:DeltaMij) into eq. (4.74), comparing with eq. (5.3), and setting for simplicity $\Delta M_{+}=\Delta M_{\bar{0}}\equiv\Delta M$, the effects from $\Delta M$ can phenomenologically be included through the substitutions

[TABLE]

We adopt this recipe in the following, setting for illustration $\Delta M=10^{-3}M$, which is parametrically in the correct range $\sim\alpha^{2}M$ and numerically in fair accordance with ref. [19] at $\lambda_{i}=0$, and also reflects the gradual increase of $\Delta M\simeq\lambda_{4,5}v^{2}/M$ with scalar self-couplings. The case $\Delta 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 $E^{\prime}\in(E^{\prime}_{\mbox{\scriptsize min}},E^{\prime}_{\mbox{\scriptsize max}})$, where E^{\prime}_{\mbox{\scriptsize min}}\equiv 2\Delta M_{\mbox{\tiny\rm{T}}}-15\alpha^{2}M and $E^{\prime}_{\mbox{\scriptsize max}}\equiv 15T$, where $\alpha\equiv(g_{1}^{2}+3g_{2}^{2})/(16\pi)$.

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_{\mbox{\scriptsize eq}}^{2}\to 2Y_{\mbox{\scriptsize eq}}(Y-Y_{\mbox{\scriptsize 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_{\mbox{\scriptsize eq}}$ has become exponentially small by the time that $z\sim 40$. In the absence of $Y_{\mbox{\scriptsize eq}}$, eq. (5.2) can be integrated into

[TABLE]

The regime z\mathop{\raise 1.29167pt\hbox{>\kern-7.5pt\raise-4.73611pt\hbox{\sim}}}40 can easily reduce the dark matter abundance by a factor $2...3$. We choose $z_{\mbox{\scriptsize 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\simeq 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 factor $\bar{S}_{1,{\mbox{\scriptsize 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 $\pi T\ll\alpha^{2}M$. The value $z_{\mbox{\scriptsize 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_{\mbox{\scriptsize phys}}=Y(z_{\mbox{\scriptsize final}})$. The energy density carried by the lightest ones today is $\rho_{\mbox{\scriptsize dm}}(T_{0})=MY_{\mbox{\scriptsize phys}}s(T_{0})$, and the energy fraction is $\Omega_{\mbox{\scriptsize dm}}(T_{0})=MY_{\mbox{\scriptsize phys}}s(T_{0})/\rho_{\mbox{\scriptsize cr}}(T_{0})$, where $\rho_{\mbox{\scriptsize cr}}$ is the current critical energy density. Inserting from ref. [58] $s(T_{0})=2\,891$/cm3 and $\rho_{\mbox{\scriptsize cr}}(T_{0})=1.0537\times 10^{-5}h^{2}\,\mbox{GeV}$/cm3 yields

[TABLE]

which can be compared with the observed value $\left.\Omega_{\mbox{\scriptsize dm}}h^{2}\right|_{\mbox{\scriptsize obs}}=0.1186(20)$ [59]. Results are plotted in fig. 6; a discussion is deferred to the first paragraph of sec. 6.

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\mathop{\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}}535\pm 9 GeV can typically be found in literature (cf. e.g. refs. [25, 36]), and that for this case we get M\mathop{\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}}519\pm 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\mathop{\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}}523\pm 5 GeV for $\Delta M=10^{-3}M$, and to M\mathop{\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}}562\pm 5 GeV for the extreme case $\Delta M/M\to 0$. For the maximal quartic couplings considered, $\lambda_{3}(2M)=\lambda_{4}(2M)=\lambda_{5}(2M)=\pi$, we obtain M\mathop{\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}}10.6\pm 0.1 TeV with free spectral functions; M\mathop{\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}}11.1\pm 0.1 TeV for $\Delta M=10^{-3}M$; and M\mathop{\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}}12.5\pm 0.1 TeV for $\Delta M/M\to 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\ll\alpha^{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 $\pi T\ll\alpha^{2}M$, there is gradual departure from kinetic and chemical equilibrium in the dark sector, whose study represents a complicated but interesting open problem.101010A nice recent investigation of non-equilibrium effects in another context can be found in ref. [60].

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, \Delta M_{\mbox{\tiny\rm{T}}}\simeq(2\lambda_{3}+\lambda_{4})T^{2}/(24M). In addition, at T\mathop{\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}}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.

Beyond the IDM, our interest lies in models including strongly interacting particles, which have attracted much recent interest in view of the substantial role that bound states could play (cf. e.g. refs. [2, 9, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78]). Having now “calibrated” the formalism of ref. [2] through a much-studied test case, we hope to address such models in the near future.

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.

Appendix Appendix A Explicit expressions for static potentials

We present here the explicit expressions for the potentials appearing in eqs. (4.17)–(4.19). The potentials contain the Debye mass parameters defined in eq. (4.16), and the charged and neutral gauge mass parameters [52]

[TABLE]

The mixing angles are defined by

[TABLE]

The neutral gauge field components are diagonalized as in eqs. (A.21)–(A.23) of ref. [2]. We define the functions

[TABLE]

Denoting furthermore $\tilde{c}\equiv\cos\tilde{\theta}$, $\tilde{s}\equiv\sin\tilde{\theta}$, $c\equiv\cos\theta$, and $s\equiv\sin\theta$, and renormalizing $r$-independent linear divergences of the potentials as mentioned below eq. (4.19), we get

[TABLE]

The potentials get considerably simplified in the short-distance limit r\ll 1/m_{\mbox{\tiny\rm{\widetilde{Z}}}}. Then their divergent $r$-dependent parts read

[TABLE]

and eqs. (4.17)–(4.19) become

[TABLE]

Defining $\alpha_{1}\equiv(3g_{2}^{2}+g_{1}^{2})/(16\pi)$, $\alpha_{2}\equiv(g_{2}^{2}-g_{1}^{2})/(16\pi)$ and $\alpha_{3,4}\equiv(g_{2}^{2}+g_{1}^{2})/(16\pi)$, the corresponding Sommerfeld factors read [10]

[TABLE]

where $X_{i}\equiv\pi\alpha_{i}/v$ and $v$ parametrizes $E^{\prime}$ from eq. (4.10) as E^{\prime}=2\Delta M_{\mbox{\tiny\rm{T}}}+Mv^{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\ll\pi T$. If m\mathop{\raise 1.29167pt\hbox{>\kern-7.5pt\raise-4.73611pt\hbox{\sim}}}\pi 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_{\mbox{\tiny\rm{E1}}}^{2}$ and $m_{\mbox{\tiny\rm{E2}}}^{2}$. 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_{\mbox{\scriptsize\sl f\,}}$; and model $W^{\pm}$, $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

[TABLE]

Here the fermionic and bosonic susceptibilities read

[TABLE]

where $n_{\mbox{\tiny\rm{F{}}}}$ and $n_{\mbox{\tiny\rm{B{}}}}$ are the Fermi and Bose distributions, respectively. We have adopted eq. (A.17) for modelling the low-temperature regime, inserting $m_{\mbox{\scriptsize\sl f\,}}\simeq(m_{t}m_{b})^{1/2}$ and m_{g}\simeq(m_{\mbox{\tiny\rm{Z}}}m_{\mbox{\tiny\rm{W}}}^{2}m_{\phi})^{1/4}, but stress that this represents a purely phenomenological recipe within the complicated temperature interval m_{b}\mathop{\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}}\pi T\mathop{\raise 1.29167pt\hbox{<\kern-7.5pt\raise-4.73611pt\hbox{\sim}}}m_{t}.