Early kinetic decoupling of dark matter: when the standard way of calculating the thermal relic density fails

TL;DR

This paper challenges the standard assumption of local thermal equilibrium during dark matter freeze-out, introducing new methods to accurately compute relic density when early kinetic decoupling occurs, significantly affecting abundance predictions.

Contribution

It presents two novel approaches to handle early kinetic decoupling in dark matter relic density calculations, improving accuracy over traditional methods.

Findings

Early kinetic decoupling can alter relic density predictions by up to an order of magnitude.

The proposed methods provide more accurate relic abundance calculations in scenarios with non-standard thermal histories.

Application to Scalar Singlet dark matter demonstrates the impact of early decoupling on abundance estimates.

Abstract

Calculating the abundance of thermally produced dark matter particles has become a standard procedure, with sophisticated methods guaranteeing a precision that matches the percent-level accuracy in the observational determination of the dark matter density. Here, we point out that one of the main assumptions in the commonly adopted formalism, namely local thermal equilibrium during the freeze-out of annihilating dark matter particles, does not have to be satisfied in general. We present two methods for how to deal with such situations, in which the kinetic decoupling of dark matter happens so early that it interferes with the chemical decoupling process: i) an approximate treatment in terms of a coupled system of differential equations for the leading momentum moments of the dark matter distribution, and ii) a full numerical solution of the Boltzmann equation in phase-space. For illustration, we apply these methods to the case of Scalar Singlet dark matter. We explicitly show that even in this simple model the prediction for the dark matter abundance can be affected by up to one order of magnitude compared to the traditional treatment.

Equations

$$E\left({\partial }_t-H\mathbf{p}\cdot {\nabla }_{\mathbf{p}}\right)f_{\chi }=C[f_{\chi }].$$

$$C_{\mathrm{ann}}$$

$$C_{\mathrm{el}}$$

$$C_{\mathrm{ann}}$$

$$C_{\mathrm{el}}\simeq\frac{m_{\chi}}{2}\gamma(T){\Bigg{[}}Tm_{\chi}\partial_{p}^{2}+\left(p+2T\frac{m_{\chi}}{p}\right)\partial_{p}+3{\Bigg{]}}f_{\chi}\,,$$

$$\gamma (T)=\frac{1}{48{\pi }^3{g}_{\chi }{m}_{\chi }^3}\int d\omega g^{\pm }{\partial }_{\omega }\left(k^4{⟨{\left|\mathcal{M}\right|}^2⟩}_t\right),$$

$${⟨{\left|\mathcal{M}\right|}^2⟩}_t\equiv \frac{1}{8k^4}{\int }_{-4k_{\mathrm{cm}}^2}^{0}dt(-t){\left|\mathcal{M}\right|}^2=16\pi m_{\chi }^2{\sigma }_T,$$

$$C_{\mathrm{el}}\simeq\frac{E}{2}\gamma(T){\Bigg{[}}TE\partial_{p}^{2}+\left(p+2T\frac{E}{p}+T\frac{p}{E}\right)\partial_{p}+3{\Bigg{]}}f_{\chi}$$

$$\frac{d{n}_{\chi }}{dt}+3H{n}_{\chi }=g_{\chi }\int \frac{d^{3}p}{(2\pi {)}^{3}E}{C}_{\mathrm{ann}}[f_{\chi }],$$

$$n_{\chi }=g_{\chi }\int d^{3}p/(2\pi {)}^3{f}_{\chi }(\mathbf{p})$$

$$f_{\chi }=A(T)f_{\chi ,\mathrm{eq}}=\frac{n_{\chi }}{n_{\chi ,\mathrm{eq}}}{f}_{\chi ,\mathrm{eq}},$$

$$\frac{d{n}_{\chi }}{dt}+3H{n}_{\chi }=⟨\sigma v⟩\left(n_{\chi ,\mathrm{eq}}^2-n_{\chi }^2\right),$$

$$⟨\sigma v⟩$$

$$x$$

$$Y$$

$$\frac{Y^{\prime }}{Y}=\frac{sY}{x\tilde{H}}⟨\sigma v⟩\left[\frac{Y_{\mathrm{eq}}^2}{Y^2}-1\right].$$

$$\tilde{g}\equiv \frac{1}{3}\frac{T}{g_{\text{eff}}^s}\frac{d{g}_{\text{eff}}^s}{dT}.$$

$${\Omega }_{\chi }{h}^2=2.755\times 10^{10}\left(\frac{m_{\chi }}{100\mathrm{GeV}}\right){\left(\frac{T_{\mathrm{CMB}}}{2.726\mathrm{K}}\right)}^3{Y}_0.$$

$$y\equiv \frac{m_{\chi }}{3s^{2/3}}⟨\frac{{\mathbf{p}}^2}{E}⟩=\frac{m_{\chi }}{3s^{2/3}}\frac{g_{\chi }}{n_{\chi }}\int \frac{d^{3}p}{(2\pi {)}^3}\frac{{\mathbf{p}}^2}{E}{f}_{\chi }(\mathbf{p}).$$

$$T_{\chi }=y{s}^{2/3}/m_{\chi }.$$

$$\frac{Y^{\prime }}{Y}$$

$$\frac{y^{\prime }}{y}$$

$$⟨p^4/E^3⟩\equiv n_{\chi }^{-1}{g}_{\chi }\int \frac{d^{3}p}{(2\pi {)}^3}\frac{{\mathbf{p}}^4}{E^3}{f}_{\chi }(\mathbf{p})$$

$$m_{\chi }{n}_{\chi }{C}_0$$

$$m_{\chi }{n}_{\chi }⟨\frac{{\mathbf{p}}^2}{E}⟩C_2$$

$$\frac{Y^{\prime }}{Y}$$

$$\frac{y^{\prime }}{y}$$

$${⟨\sigma v⟩}_2$$

$${⟨\sigma v⟩}_{\mathrm{neq}}$$

$${⟨\sigma v⟩}_{2,\mathrm{neq}}$$

$$⟨p^4/E^3{⟩}_{\mathrm{neq}}$$

$$f_{\chi}=\frac{n_{\chi}(T)}{n_{\chi,\text{eq}}(T_{\chi})}\exp\left(-\frac{E}{T_{\chi}}\right)\bigg{|}_{T_{\chi}=ys^{2/3}/m_{\chi}}\,,$$

$$T_{\chi }\left[\frac{y_{\mathrm{eq}}}{y}-1\right]=T-T_{\chi }$$

Full text

I Introduction

The leading hypothesis for the cosmological dark matter (DM) Hinshaw et al. (2013); Ade et al. (2016) is a new type of elementary particle Bertone et al. (2005). One of the most attractive options to explain the present abundance of these particles consists in the possibility that they have been thermally produced in the early universe. This is particularly interesting for the scenario originally studied by Lee and Weinberg Lee and Weinberg (1977), as well as others Hut (1977); Sato and Kobayashi (1977); Dicus et al. (1977); Wolfram (1979), in which non-relativistic DM particles initially are kept in thermal equilibrium with the heat bath through frequent annihilation and creation processes with standard model (SM) particles. Once the interaction rate starts to fall behind the expansion rate of the universe, the DM number density begins to ‘freeze out’ and remains covariantly conserved. For weakly interacting massive particles (WIMPs), elementary particles with masses and interaction strengths at the electroweak scale, this scenario automatically leads to a relic abundance in rough agreement with the observed DM density — a fact sometimes referred to as the WIMP miracle.

The by now standard treatment Gondolo and Gelmini (1991); Edsjo and Gondolo (1997) of calculating the resulting DM abundance in these scenarios implements an efficient and highly accurate method of solving the Boltzmann equation for a given (effective) invariant DM annihilation rate. This approach fully captures, in particular, the three famous exceptions to the original relic density calculations pointed out in a seminal paper by Griest and Seckel Griest and Seckel (1991), namely co-annihilations, threshold effects and resonances. The main assumption entering this formalism is that, during the freeze-out process, DM is still kept in local thermal equilibrium with the heat bath by frequent scattering processes with relativistic SM particles. For many WIMP candidates, this is indeed satisfied to a high accuracy and kinetic decoupling typically only happens much later than the chemical decoupling Bringmann (2009).

Here we point out that exceptions to this standard lore do exist, even in very simple scenarios, where kinetic decoupling happens so early that it cannot be neglected during the freeze-out process. We develop both semi-analytical and fully numerical methods to solve the Boltzmann equation and to compute the DM relic abundance in these cases. Technically, one of the challenges that had to be overcome for obtaining sufficiently accurate results was to extend the highly non-relativistic Boltzmann equation, as discussed previously in the literature, to the semi-relativistic regime. Numerically, we also succeeded to resolve the evolution of the full phase-space distribution accurately enough to test, for the first time, the underlying assumptions for the standard way of calculating the relic density of WIMPs or other self-annihilating DM candidates (for a recent example where the relic density is instead set by inelastic scattering, rather than self-annihilation, see Ref. D’Agnolo et al. (2017)). We illustrate our general results by a detailed discussion of the Scalar Singlet model Silveira and Zee (1985); McDonald (1994); Burgess et al. (2001), for which we find a DM relic density that differs by up to an order of magnitude from the standard treatment.

This article is organized as follows. In Section II, we start with a general description of the underlying Boltzmann equation that governs the DM phase-space evolution. We then briefly review the standard treatment of solving for the DM number density (Section II.1), extend this by deriving a coupled system of evolution equations for the number density and the velocity dispersion (Section II.2), and finally introduce our framework for a fully numerical solution (Section II.3). Section III is devoted to a thorough application of these methods to the Scalar Singlet model. We comment on our results in Section IV, and discuss potential other areas of application, before we conclude in Section V. In two Appendices we discuss in detail the evolution of the Singlet DM phase-space density for selected parameter points (App. A) and comment on the semi-relativistic form of the scattering operator in the Boltzmann equation (App. B).

II Thermal production of dark matter

Let us denote the DM particle by $\chi$, and its phase-space density by $f_{\chi}(t,\mathbf{p})$. The evolution of $f_{\chi}$ is governed by the Boltzmann equation which, in an expanding Friedmann-Robertson-Walker universe, is given by Kolb and Turner (1990); Bringmann and Hofmann (2007)

[TABLE]

Here, $H=\dot{a}/a$ is the Hubble parameter, $a$ the scale factor, and the collision term $C[f_{\chi}]$ contains all interactions between DM and SM particles $f$. For WIMPs, we are to leading order interested in two-body processes for DM annihilation and elastic scattering, $C=C_{\rm ann}+C_{\rm el}$, where

[TABLE]

and

[TABLE]

In the above expressions, ${\left|\mathcal{M}\right|}^{2}$ refers to the respective squared amplitude, summed over all spin and other internal degrees of freedom, as well as all SM particles $f$. We assume the SM particles to be in thermal equilibrium, such that their phase-space distribution is given by $g^{\pm}(\omega)=1/\left[\exp(\omega/T)\pm 1\right]$. Note that we have neglected Bose enhancement and Pauli blocking factors for $f_{\chi}$ here, as we assume DM to be nonrelativistic; momentum conservation then implies that, in $C_{\mathrm{ann}}$, we can also neglect these factors for the SM particles.

Assuming $CP$ invariance, and using the fact that in thermal equilibrium annihilation and creation processes should happen with the same frequency, the annihilation term given by Eq. (II) can be further simplified to Gondolo and Gelmini (1991)

[TABLE]

where $v=v_{\rm M\o l}\equiv({E\tilde{E}})^{-1}[{(p\cdot\tilde{p})^{2}-m_{\chi}^{4}}]^{1/2}$ is the Møller velocity, which in the rest frame of one of the DM particles coincides with the lab velocity $v_{\rm lab}=[s(s-4m_{\chi}^{2})]^{1/2}/{(s-2m_{\chi}^{2})}$.

The scattering term, on the other hand, is in general considerably more difficult to manage. Analytic expressions have, however, been obtained in the highly non-relativistic limit of the DM particles, and assuming that the momentum transfer in the scattering process is much smaller than the DM mass Bringmann and Hofmann (2007); Bringmann (2009); Kasahara (2009); Gondolo et al. (2012); Binder et al. (2016); Bringmann et al. (2016):

[TABLE]

where the momentum exchange rate is given by

[TABLE]

with

[TABLE]

and ${k}_{\mathrm{cm}}^{2}\!=\!\left(s-(m_{\chi}-m_{f})^{2}\right)\left(s-(m_{\chi}+m_{f})^{2}\right)/(4s)$ evaluated at ${s=m_{\chi}^{2}+2\omega m_{\chi}+m_{f}^{2}}$. Here, $\sigma_{T}=\int d\Omega(1-\cos\theta)d\sigma/d\Omega$ is the standard transfer cross section for elastic scattering. In Appendix B, we discuss how the scattering term is expected to change in the semi-relativistic case, i.e. when the assumption of highly non-relativistic DM is slightly relaxed. For reference, we will in the following use

[TABLE]

when explicitly addressing this regime.

II.1 The standard treatment

In order to calculate the DM relic abundance, we can integrate the Boltzmann Eq. (1) over $\mathbf{p}$. This results in

[TABLE]

which has to be solved for the DM number density

[TABLE]

(note that $C_{\rm el}$ vanishes once it is integrated over). In order to evaluate the r.h.s. of this equation, the usual assumption Gondolo and Gelmini (1991) is that during chemical freeze-out one can make the following ansatz for the DM distribution:

[TABLE]

where $A(T)=1$ in full equilibrium, i.e. before chemical freeze-out. This is motivated by the fact DM-SM scattering typically proceeds at a much faster rate than DM-DM annihilation, because the number density of relativistic SM particles is not Boltzmann suppressed like that of the non-relativistic DM particles. In that case, DM particles are kept in local thermal equilibrium even when the annihilation rate starts to fall behind the Hubble expansion and chemical equilibrium can no longer be maintained.

Approximating furthermore $f_{\chi,{\rm eq}}(E)\simeq\exp(-E/T)$, i.e. neglecting the impact of quantum statistics for non-relativistic particles, five of the six integrals in Eq. (9) can be performed analytically. This by now standard treatment, as established by Gondolo & Gelmini Gondolo and Gelmini (1991), results in the often-quoted expression

[TABLE]

where $n_{\chi,{\rm eq}}=g_{\chi}m_{\chi}^{2}TK_{2}(m_{\chi}/T)/(2\pi^{2})$ and

[TABLE]

Here, $K_{i}$ are the modified Bessel functions of order $i$, and we have introduced $\tilde{s}\equiv s/(4m_{\chi}^{2})$. While there are various ways to state the final result for $\left\langle\sigma v\right\rangle$, the form given above stresses that physically one should indeed think of this quantity as a thermal average of $\sigma v_{\rm lab}$ rather than any other combination of cross section and velocity (in the sense that we strictly have $\left\langle\sigma v\right\rangle=\sigma v_{\rm lab}$ for $\sigma v_{\rm lab}=$ const; for e.g. $\sigma v_{\rm CMS}=$ const, on the other hand, with $v_{\rm CMS}=2\sqrt{1-4m_{\chi}^{2}/s}$ being the relative velocity in the CMS frame, we instead have $\left\langle\sigma v\right\rangle\to\sigma v_{\rm CMS}$ only in the limit $T\to 0$).

By introducing dimensionless variables

[TABLE]

and assuming entropy conservation, finally, the above Boltzmann equation for the number density, Eq. (12), can be brought into an alternative form that is particularly suitable for numerical integration:

[TABLE]

Here, $s=(2\pi^{2}/45)g^{s}_{\rm{eff}}T^{3}$ denotes the entropy density, ${}^{\prime}\equiv d/dx$ and $\tilde{H}\equiv H/\left[1+\tilde{g}(x)\right]$ where

[TABLE]

The value of $Y$ today, $Y_{0}\equiv Y(x\to\infty)$, can then be related to the observed DM abundance by Gondolo and Gelmini (1991)

[TABLE]

We note that Eq. (17) is the basis for the implementation of relic density calculations in all major numerical codes Belanger et al. (2002); Gondolo et al. (2004); Belanger et al. (2007); Arbey and Mahmoudi (2010); Belanger et al. (2014); Bringmann et al. (2018, 2017).

II.2 Coupled Boltzmann equations

The main assumption that enters the standard treatment reviewed above is contained in Eq. (11), i.e. the requirement that during chemical freeze-out, or in fact during any period when the comoving DM density changes, local thermal equilibrium with the heat bath is maintained. If that assumption is not justified, one has in principle to solve the full Boltzmann equation in phase space, Eq. (1), numerically (see next subsection). As first pointed out in Ref. van den Aarssen et al. (2012), however, it sometimes suffices to take into account the second moment of Eq. (1), instead of only the zeroth moment as in the previous subsection. This leads to a relatively simple coupled system of differential equations that generalizes Eq. (17).

The starting point is to define, in analogy to $Y$ for the zeroth moment of $f_{\chi}$, a dimensionless version of the second moment of $f_{\chi}$:

[TABLE]

For a thermal distribution, the DM particles thus have a temperature

[TABLE]

We note that for non-thermal distributions we could still view this last equation as an alternative definition of the DM ‘temperature’, or velocity dispersion, in terms of the second moment of $f_{\chi}$ as introduced above. This allows, e.g., a convenient characterization of kinetic decoupling as the time when $T_{\chi}$ no longer equals $T$ but instead starts to approach the asymptotic scaling of $T_{\chi}=T_{\rm kd}(a/a_{\rm eq})^{-2}$ for highly non-relativistic DM Bringmann and Hofmann (2007); Bringmann (2009).

Integrating Eq. (1) over $g_{\chi}\int d^{3}p/(2\pi)^{3}/E$ and $g_{\chi}\int d^{3}p/(2\pi)^{3}\mathbf{p}^{2}/E^{2}$, respectively, we find

[TABLE]

where

[TABLE]

and we introduced the moments of the collision term as

[TABLE]

Plugging in $C=C_{\rm ann}+C_{\rm el}$ as provided in Eqs. (4,5), finally, we arrive at a coupled set of equations that constitutes one of our main results:111This extends the results presented in van den Aarssen et al. (2012). Compared to that reference, we have kept terms proportional to $Y_{\rm eq}$ (see also Duch and Grzadkowski (2017)) and adopted a fully relativistic temperature definition in Eqs. (20,21). The latter indeed turns out to be important outside the highly non-relativistic regime and is the origin of the last term in Eq. (II.2), as well as the corrected form of $\left\langle\sigma v\right\rangle_{2}$ – which now (unlike in its original form) can be seen as a proper thermal average in the sense that a constant $\sigma v_{\rm lab}$ leads to $\left\langle\sigma v\right\rangle_{2}=\sigma v_{\rm lab}$ for all values of $T$ (i.e. not only for $T\to 0$).

We note that both $\left\langle p^{4}/E^{3}\right\rangle$ and the integral over $\epsilon_{+}$ can be expressed in terms of a series of Bessel functions when expanding $E$ in the denominator around $E=m$. Since this series does not converge very fast for the relatively small values of $x$ that we will be interested in here, however, we do not display these series.

[TABLE]

Here, in addition to $\langle\sigma v\rangle$ in Eq. (II.1), we also introduced another, temperature-weighted thermal average:

[TABLE]

where we have used $\epsilon_{+}\equiv(E+\tilde{E})/\sqrt{s}$. The ‘out-of equilibrium average’ $\left\langle\sigma v\right\rangle_{2,{\rm neq}}$ is defined as in Eq. (30), but for arbitrary $n_{\chi}$, $f_{\chi}(\mathbf{p})$ – and hence also $1/T\to 1/T_{\chi}$ in the normalization; the last equality, Eq. (30), thus does not hold in this case. Correspondingly, $\left\langle\sigma v\right\rangle_{\rm neq}$ is defined in analogy to Eq. (II.1), but equals in general not the expression given in Eq. (II.1).

Two comments about this central result are in order. The first comment, more important from a practical point of view, is that the set of equations (22, 23) includes higher moments of $f_{\chi}$, and hence does not close w.r.t. the variables $Y$ and $y$. Concretely, we need additional input to determine the quantities $\left\langle\sigma v\right\rangle_{{\rm neq}}$, $\left\langle\sigma v\right\rangle_{2,{\rm neq}}$ and $\langle p^{4}/E^{3}\rangle_{\rm neq}$ in Eqs. (27, II.2) in terms of only $y$ and $Y$. We will make the following ansatz for these quantities:

[TABLE]

These expressions would, in particular, result from a DM phase-space distribution of the form

[TABLE]

which describes a situation in which the DM particles follow a Maxwellian velocity distribution with a temperature different from that of the heat bath (as expected, e.g., if the DM particles exhibit significant self-scattering Feng et al. (2010a); Buckley and Fox (2010); Feng et al. (2010b); van den Aarssen et al. (2012)). We emphasize, however, that from the point of view of solving the coupled set of equations (27, II.2), there is no need to make such a relatively strong assumption about $f_{\chi}(\mathbf{p})$: any form of $f_{\chi}$ that leads to (very) similar results for the quantities given in Eqs. (31) – (33) will also lead to (very) similar results for $Y(x)$ and $y(x)$. In other words, we expect our coupled system of Boltzmann equations to agree with the full numerical solution discussed in the next section – concerning the evolution of $Y$ and $y$ – if and only if the ansatz in Eqs. (31,32,33) coincides with the corresponding averages numerically determined from the ‘true’ phase-space distribution. As we will see later, this is indeed very often the case.

The second comment concerns the first term on the r.h.s. of Eq. (II.2), which is proportional to the second moment of the elastic scattering term given in Eq. (5). As that latter expression is valid only to lowest order in $p^{2}/E^{2}\sim p^{2}/m_{\chi}^{2}\sim 1/x$, we had for consistency also to neglect any higher-order corrections in these quantities to the elastic scattering part of $C_{2}$ when deriving our final result. As discussed in Appendix B, in fact, there is no simple way of determining the next-to-leading order corrections to $C_{\rm el}$. If we use our default semi-relativistic scattering term given in Eq. (8), however, including the resulting corrections from sub-leading orders corresponds to replacing in Eq. (II.2)

[TABLE]

By construction, see Appendix B, this operator must still be an attractor to the equilibrium solution, and hence be proportional to (some power of) $T-T_{\chi}$; for the ansatz of Eq. (34), e.g., this can easily be verified directly. In practice, this replacement has very little impact on the evolution of $Y$ and $y$, even at times as early as $x\sim 10$. We can think of the resulting small differences as a measure of the intrinsic uncertainty associated to our treatment of the scattering term.

II.3 The full phase-space density evolution

We now turn to solve the Boltzmann Eq. (1) at the full phase-space density level. This is numerically more challenging, but allows to assess the validity of the assumptions in previous sections and to track deviations (as we will see can occur) from the standard Maxwell Boltzmann velocity distribution. To achieve this, we start by re-expressing Eq. (1) in the two dimensionless coordinates

[TABLE]

where the monotonic temperature $T(t)$ replaces as before the time parameter $t$ via our $x(T)$, and $q$ is now the ‘momentum’ coordinate that depends on both $t$ and $p$. In these variables, we can rewrite the Liouville operator on the l.h.s. of Eq. (1) as

[TABLE]

Here, we used the fact that the system is isotropic and assumed, as in the previous sections, that entropy is conserved. With the collision terms given in Eqs. (4) and (8), the Boltzmann equation for $f_{\chi}$ now becomes

[TABLE]

where $x_{q}\equiv\sqrt{x^{2}+q^{2}}$ and $\theta$ is the angle between $\bf{q}$ and $\bf{\tilde{q}}$.

The benefits of this rewriting is two-fold. First, the interpretation of the Boltzmann equation becomes very transparent, in the sense that this “comoving” phase space density $f_{\chi}(x,q)$ clearly stays unaltered for $\tilde{g}(x)=0$ and vanishing annihilation and scattering rates (being proportional to $\sigma_{\bar{\chi}\chi\rightarrow\bar{f}f}$ and $\gamma$, respectively). The new coordinates thus absorb how momentum and DM density change exclusively due to the Hubble expansion. (For non-vanishing $\tilde{g}$, these quantities continue to scale in the same way with the scale factor $a$, but taking into account that $a\propto{g_{\rm eff}^{s}}^{-1/3}T^{-1}$). Second, the use of a comoving momentum $q\equiv p/T$ significantly helps numerical calculations that extend over a large range in $x=m_{\chi}/T$. In fact, $f_{\chi}(x,q)$ is expected to stay unchanged in shape both in the early semi-relativistic and kinetically coupled regime, where $f_{\chi}\sim e^{-p/T_{\chi}}=e^{-q}$ given that $T_{\chi}=T$, as well as in the late non-relativistic kinetically decoupled regime, where $f_{\chi}\sim e^{-p^{2}/(2mT_{\chi})}\propto e^{-q^{2}/(2m)}$ given that $T_{\chi}\propto T^{2}$ in this case — at least as long as $\tilde{g}=0$ and the DM phase-space distribution remains close to Maxwellian as in Eq. (34).

Let us stress that here, unlike for our discussion in the previous subsection, it is indeed mandatory to use the semi-relativistic form of Eq. (8) for the scattering operator when discussing the evolution of the phase-space density, in the sense that it must drive the distribution function $f_{\chi}(q)$ towards the fully relativistic form $\propto e^{-E/T}$ (and not as Eq. (5) to the non-relativistic approximation $\propto e^{-\frac{q^{2}}{2m}/{T}}$). The importance of this can be seen by comparing the second and the third line of Eq. (37). The term in the second line will always drive DM annihilation to occur unless an equilibrium distribution $f_{\rm eq}$ is reached. The term in the third line determines towards which equilibrium shape the scattering operator will drive the DM distribution $f_{\chi}(q)$. If the scattering attractor distribution would not match the $f_{\rm eq}(q)$ of the second line, then scattering could artificially drive annihilation to occur. For more discussions of the semi-relativistic aspects of the scattering term, see Appendix B.

We then use a technique that discretizes the unbounded momentum variable $q$ into a finite number of $q_{i}$ with $i\in\{1,2,\ldots,N\}$. This allows to rewrite our integro partial differential equation into a set of $N$ coupled ordinary differential equations (ODEs):

[TABLE]

where $f_{i}\equiv f_{\chi}(x,q_{i})$, and the derivatives $\partial_{q}f_{i}$ and $\partial^{2}_{q}f_{i}$ are determined numerically by finite differentials using several neighboring points to $f_{i}$. $\langle v_{\rm M\o l}\sigma_{\bar{\chi}\chi\rightarrow\bar{f}f}\rangle_{i,j}^{\theta}$ is the velocity-weighted cross section averaged over $\theta$ (which is evaluated analytically or numerically) as a function of $q_{i}$ and $\tilde{q}_{j}$, and $\Delta\tilde{q}_{j}\equiv\tilde{q}_{j+1}-\tilde{q}_{j}$. Finally, the DM number density in Eq. (10) is determined by trapezoidal integration.

Numerous numerical tests have been performed to ensure stability of our solutions to the ODEs of Eq. (38) and that imposed conditions on the now emerged boundary points (at $q_{1}$ and $q_{N}$) are physically sound. It turns out that very small stepsizes over a large range in $q$ are required for solving these stiff ODEs. We typically used the range $q_{1}=10^{-6}$ to $q_{N}=50$ with about thousand steps in between, and set the two last terms of Eq. (38) to zero at $q_{N}$ while using forward derivatives to evaluate them at $q_{1}$. By the use of the ODE15s code in MatLab, and by analytically deriving internally required Jacobians, we are able to efficiently calculate the full phase-space evolution for the freeze-out after optimizing numerical settings. On the time scale of a few minutes we can derive the relic abundance for a given DM model. The code is general enough to be adapted to any standard single WIMP setup.

III Scalar Singlet Dark Matter

The simplest example of a renormalizable model providing a WIMP DM candidate is the Scalar Singlet model Silveira and Zee (1985); McDonald (1994); Burgess et al. (2001), originally proposed as DM made of ‘scalar phantoms’ by Silveira and Zee Silveira and Zee (1985). In this model, the only addition to the standard model is a real gauge-singlet scalar field $S$ which is stabilized by a $\mathbb{Z}_{2}$ symmetry and never obtains a non-vanishing vacuum expectation value. The simplicity of the model has in itself triggered considerable interest Yaguna (2009); Profumo et al. (2010); Arina and Tytgat (2011); Mambrini (2011); Lerner and McDonald (2009); Herranen et al. (2015); Kahlhoefer and McDonald (2015); Profumo et al. (2007); Barger et al. (2009); Cline and Kainulainen (2013), with a further boost of attention after the discovery of the Higgs boson Djouadi et al. (2012); Cheung et al. (2012); Endo and Takaesu (2015); Djouadi et al. (2013); Cline et al. (2013); Urbano and Xue (2015); He and Tandean (2016); Escudero et al. (2016); Goudelis et al. (2009); Craig et al. (2016); Han et al. (2016); Ko and Yokoya (2016); Beniwal et al. (2016); Cuoco et al. (2016). Recently, the GAMBIT Athron et al. (2017a) collaboration presented the so far most comprehensive study of this model by performing a global fit taking into account experimental constraints from both direct, indirect and accelerator searches for DM Athron et al. (2017b).

Interestingly, the resulting parameter region with the highest profile likelihood in this global fit is the one where the Scalar Singlet mass $m_{S}$ is about half that of the SM Higgs mass, $m_{h}$, and where the DM abundance today is set by the resonant annihilation of two DM particles through an almost on-shell Higgs boson. As we will see, it is exactly in this parameter region that the standard way of calculating the relic density, as implemented in all previous studies of this model, fails because kinetic decoupling happens so early that it essentially coincides with chemical decoupling. Instead, the formalism introduced in the previous section provides a reliable calculation of the relic abundance of Scalar Singlet DM.

III.1 Model setup

The model symmetries, along with the requirement of renormalizability, uniquely determine the form of the Lagrangian to be

[TABLE]

where $H$ is the SM Higgs doublet. The $S$ boson mass receives contributions from both the bare mass term, $\mu_{S}$, and from electroweak symmetry breaking, leading to $m_{S}=\sqrt{\mu_{S}^{2}+\frac{1}{2}{\lambda_{S}v_{0}^{2}}}$, where $v_{0}=246.2$ GeV is the Higgs vacuum expectation value. We adopt the Higgs mass and the total width from decay to SM particles to be $m_{h}=125.09$ GeV Patrignani et al. (2016) and $\Gamma_{h,SM}=4.042$ MeV Dittmaier et al. (2011). For the moment, we neglect the quartic self-coupling $\lambda_{SS}$, but will later comment on its potential (minor) impact on relic density calculations.

The annihilation cross section of DM pairs to SM particles, apart from $hh$ final states, is given by Cline et al. (2013)

[TABLE]

where $\Gamma_{h\to\rm{SM}}(\sqrt{s})$ is the partial decay width of a Standard-Model Higgs boson of mass $\sqrt{s}$, and

[TABLE]

The total Higgs width $\Gamma_{h}$ in the above propagator, but not elsewhere, includes not only all SM channels but also the $h\to SS$ channel if it is open. For $\Gamma_{h\to\rm{SM}}(\sqrt{s})$, as in Cline et al. (2013), we use tabulated values for $\sqrt{s}<300$ GeV from Dittmaier et al. (2011) and analytic expressions at higher $\sqrt{s}$. Note however that the latter high $\sqrt{s}$ region has no impact on the relic density in the studied Scalar Singlet mass range. Likewise, the channel $SS\to hh$ lies outside our kinematic region of interest.

For the elastic scattering processes, we take into account DM scattering with all SM fermions. Being mediated only by a Higgs in the $t$-channel, the corresponding squared amplitude takes a particularly simple form,

[TABLE]

where $m_{f}$ is the mass of the SM fermion and the color factor is $N_{f}=3$ for quarks and $N_{f}=1$ for leptons. Averaging over the transferred momentum, as in Eq. (7), we thus find

[TABLE]

Note that the sum here runs over all relevant fermions and antifermions separately.

The hierarchical Yukawa structure of the Higgs couplings leads to the scattering rate being dominated by the interactions with the heaviest fermions that for a given temperature are still sufficiently abundant in the plasma. In the range of DM masses $m_{S}$ that we are interested in, freeze-out happens around $T\sim\mathcal{O}(1\,\rm{GeV})$, which is not far from the temperature of the QCD phase transition. Consequently, the details of this transition and the SM plasma can have a significant impact on the scattering rate; a study which goes beyond the scope of this work. Therefore, we follow the literature and adopt two extreme scenarios that can be thought of as bracketing the actual size of the scattering term:

A

all quarks are free and present in the plasma down to temperatures of $T_{c}=154$ MeV (largest scattering scenario, as adopted in Gondolo et al. (2012))

B

only light quarks ($u$, $d$, $s$) contribute to the scattering, and only for temperatures above $4T_{c}\sim 600$ MeV, below which hadronization effects start to become sizeable Boyanovsky et al. (2006) (smallest scattering scenario, as adopted in Bringmann (2009)).

Finally, we adopt the recent results from Drees et al. Drees et al. (2015) for the effective number of relativistic degrees of freedom $g_{\rm eff}(T)$ that enter the calculation of the Hubble rate during radiation domination, $H=\sqrt{4\pi^{3}g_{\rm eff}/45}\,T^{2}/m_{\rm Pl}$, as well as the entropy degrees of freedom entering for example in the calculation of $\tilde{g}(T)$ as defined in Eq. (18).

III.2 Relic density of scalar singlet dark matter

Let us first compute the relic density following the standard treatment adopted in the literature. To this end, we numerically solve Eq. (17) for a given set of parameters $(m_{S},\lambda_{S})$ and determine the resulting asymptotic value of $Y_{0}$. The blue dashed line in Fig. 1 shows the contour in this plane that results in $Y_{0}$ corresponding to a relic density of $\Omega h^{2}=0.1188$, c.f. Eq. (19). We restrict our discussion to values of $m_{S}$ in the kinematic range where $\langle\sigma v\rangle$ is enhanced due to the Higgs propagator given in Eq. (41), and the coupling $\lambda_{S}$ that results in the correct relic density is hence correspondingly decreased. This curve agrees with the corresponding result obtained in Ref. Cline et al. (2013).

For comparison, we show in the same figure the required value of $\lambda_{S}$ that results when instead solving the coupled system of Boltzmann equations (27) and (II.2), or when numerically solving the full Boltzmann equation as described in Section II.3 . Here, the solid (dashed) line shows the situation for the ‘B’ (‘A’) scenario for scatterings on quarks. Outside the resonance region, the coupled Boltzmann equations lead to identical results compared to the standard approach, indicating that kinetic decoupling indeed happens much later than chemical decoupling and that the assumption of local thermal equilibrium during chemical freeze-out thus is satisfied. For DM masses inside the resonance region, on the other hand, we can see that the two methods can give significantly different results, implying that this assumption must be violated. For the same reason, a smaller scattering rate (as in scenario ‘B’) leads to an even larger deviation from the standard scenario than the maximal scattering rate adopted in scenario ‘A’.

This interpretation is explicitly confirmed in Fig. 2, where we plot the temperatures at which the DM number density and temperature start to deviate from the equilibrium values: in the parameter range that we focus on here, kinetic decoupling happens indeed very close to chemical decoupling. The reason for this very early kinetic decoupling is straight-forward to understand as the result of a strongly suppressed momentum transfer rate $\gamma(T)$, compared to the annihilation rate, due to two independent effects: i) the small coupling $\lambda_{S}$ needed to satisfy the relic density requirement, without a corresponding resonant enhancement of $\gamma(T)$, and ii) the scattering rate being proportional to the Yukawa coupling squared, which favours scattering with Boltzmann-suppressed heavy fermions. We note that the latter point also explains the relatively large difference between the two extreme quark scattering scenarios used here for illustration (in scenario ‘B’, the largest Yukawa couplings do not contribute to the scattering).

In order to emphasize the importance of our improved treatment of the decoupling history, we plot in Fig. 3 also the ratio of the resulting relic density to that of the standard approach (for parameter values satisfying the relic density constraint for the latter, i.e. corresponding to the blue dashed curve in Fig. 1). Let us stress that, compared to the observational uncertainty in this quantity of about 1 %, these corrections are by no means small even in the minimal scattering scenario ‘A’. In the same figure, we also compare our result for the coupled system of Boltzmann equations (27) and (II.2) to the full numerical solution of the Boltzmann equation in phase space, as described in Section II.3 (black dots). Before getting back to these results, let us briefly comment on the green dashed line in Fig. 3, which implements the highly non-relativistic scattering term $C_{\text{el}}$ of Eq. (5), and hence not the replacement (35) in Eq. (II.2) which we otherwise adopt as our default. Clearly, the impact of this choice is very limited for this approach. We note that the quantitative importance of the relativistic correction term proportional to $\langle p^{4}/E^{3}\rangle$ in Eq. (II.2) lies in the same ballpark, affecting the relic density by at most \sim$$10\% in the region very close to the resonance (and below the percent-level elsewhere).

In Appendix A we discuss in depth the time evolution of both the coupled Boltzmann equations and the full phase-space density in the resonance region. Let us here just mention that the characteristic features of the curves displayed in Figs. 2 and 3 can indeed all more or less directly be understood in terms of the highly enhanced annihilation rate in a relatively narrow kinematic region around the resonance, $\sqrt{s}\sim m_{h}\pm\Gamma_{h}$. As the full numerical solution reveals, furthermore, the shape of $f_{\chi}(p)$ can in some cases be quite different from the Maxwell-Boltzmann form (34) that is consistent with the coupled system of Boltzmann equations (27) and (II.2). Whether this has a noticeable impact on the resulting relic density (like for $m_{S}\sim 57$ GeV) or not (like for $m_{S}\sim m_{h}/2$) again mostly depends on whether or not the shape is affected for those momenta that can combine to $\sqrt{s}\sim m_{h}$ during chemical freeze-out.

For illustration, we pick a DM mass of $m_{S}=57$ GeV and show in Fig. 4 the full phase-space distribution for a few selected values of $x$ (left panel) as well as the relevant evolution of $Y$ and $y$ (right panel). For models with DM masses in this range, the relatively large difference between full solution and coupled equations (as visible in Fig. 3) can mostly be understood in terms of the dip in the ratio of DM phase-space distributions at intermediate values of $q=p/T$ that starts to develop for $x\gtrsim 20$. Concretely, the fact that the actual distribution for those momenta is slightly suppressed compared to a distribution fully characterized only by its second moment, as in Eq. (34), causes the DM particles to annihilate less efficiently, $\langle\sigma v\rangle_{\mathrm{neq}}<\langle\sigma v\rangle$, because this is the momentum range probed by the resonance for these $x$ values. This in turn leads to the DM particles falling out of chemical equilibrium earlier, and hence a larger asymptotic value of $Y$. The reason for this momentum suppression to develop in the first place is also to be found in the particularly efficient annihilation close to the resonance, which leads to a depletion of DM particles with corresponding momenta because the scattering rate is no longer sufficiently large to redistribute the phase-space distribution to a thermal shape. We note that the bulk part of this effect is actually well captured by the coupled Boltzmann system, c.f. the dashed vs. solid lines in the right panel of Fig. 4. For further details, we refer again to Appendix A.

IV Discussion

From the above discussion, we have learned that very early kinetic decoupling is not just a theoretical possibility. It can appear in simple WIMP models, like the Scalar Singlet case, and affect the DM relic density in a significant way. We note that the size of the latter effect is, as expected, directly related to the size of the momentum exchange rate and hence to just how early kinetic decoupling happens compared to chemical decoupling. Let us stress that, from a general point of view, this is a much more important message connected to our choice of considering two scattering scenarios than the question of which of those scenarios is more realistic for the specific model we have studied here.

We have also seen that the coupled system of Boltzmann equations (27) and (II.2) provides a qualitatively very good description for the resulting DM abundance, see in particular Fig. 1, even though for high-precision results it seems mandatory to actually solve the full Boltzmann equation in phase space. As discussed in Appendix A, differences can arise when the true phase-space distribution is not of the Maxwellian form assumed in Eq. (34) – though the two methods can actually still give almost identical results for the relic abundance even when the two distribution differ vastly. The question of under which conditions the coupled system of equations provides an accurate description of the relic density is thus a somewhat subtle one, and requires a careful discussion of the velocity dependence of the annihilation term in the Boltzmann equation.

An exception to this general complication is a DM self-interaction rate large enough to force the DM distribution into the form given by Eq. (34) Feng et al. (2010a); Buckley and Fox (2010); Feng et al. (2010b); van den Aarssen et al. (2012) and hence render the coupled system of Boltzmann equations (27) and (II.2) exactly correct (up to, as discussed, corrections due to quantum statistics). Sizeable self-scattering rates can for example arise due to corresponding contact interactions, like the quartic coupling $\lambda_{SS}$ in the Scalar Singlet case, or by adding light mediators that couple to the DM particle (which was indeed the first time such a coupled system of Boltzmann equations was considered van den Aarssen et al. (2012), albeit in a different context). For the case of resonant annihilation, furthermore, the same resonance also mediates an enhanced self-interaction. For future work, it would hence be worthwhile to extend our numerical framework to even include those DM self-interaction processes. For the Scalar Singlet case, in particular, we expect that adding the process $SS\to h^{*}\to SS$ would bring all numerical results for the full Boltzmann equation – e.g. those shown in Fig. 3 – even closer to those resulting from the coupled system of Boltzmann equations.

Let us finally stress that both the coupled Boltzmann equations and the numerical setup that we have described here are very general, and can be used to consistently study early kinetic decoupling for a much larger range of models than the Scalar Singlet case. Obvious applications are other scenarios where resonant annihilation and/or annihilation to heavy final states is important in setting the relic abundance, see also Ref. Duch and Grzadkowski (2017). Further examples where the ratio of the scattering rate to the annihilation rate can be smaller than usual, hence potentially leading to early kinetic decoupling, include Sommerfeld-enhanced annihilation Dent et al. (2010); Zavala et al. (2010); Feng et al. (2010b); van den Aarssen et al. (2012) (if the light mediators are not abundant enough to take part in the scattering process) and annihilation to DM bound states Feng et al. (2009); von Harling and Petraki (2014). Quite in general, our methods provide a powerful means to check whether the DM particles are indeed in local thermal equilibrium with the heat bath around the time when their abundance freezes out – which is the usual assumption, though rarely explicitly tested, not only in WIMP-like scenarios but also when so-called semi-annihilations D’Eramo and Thaler (2010) are important in setting the relic density, when computing the relic abundance for modified expansion histories D’Eramo et al. (2017); Redmond and Erickcek (2017), or in scenarios that go beyond simple $2\to 2$ annihilation processes Hochberg et al. (2014, 2015); Kuflik et al. (2016).

V Conclusions

The standard way of calculating the thermal relic density of self-annihilating DM particles rests on the assumption of local thermal equilibrium during freeze-out, and that hence kinetic decoupling occurs much later than chemical decoupling. Here, we demonstrated for the first time that departure from kinetic equilibrium can instead happen much earlier, even simultaneously with the departure from chemical equilibrium.

By introducing a coupled system of equations for the DM number density and its ‘temperature’, or rather velocity dispersion, we improved the standard way of calculating the relic density in such cases. For an even higher accuracy in predicting the DM abundance, we also found a way of solving the full Boltzmann equation numerically. The latter approach has the additional advantage of obtaining the full phase-space distribution, rather than only the number density, which in particular allows to test in detail the assumption of a Maxwellian velocity distribution adopted in the standard approach. A numerical solver for the coupled system of Boltzmann equations, Eqs. (27, II.2), will be available in an upcoming version of DarkSUSY Bringmann et al. (2018) and our implemented solver for the full Boltzmann equation at the phase-space level, Eq. (38), will be released separately.222Please contact any of the authors if you need these numerical routines prior to their public release.

Applied to the simplest renormalizable WIMP model – the Scalar Singlet, extensively discussed in the literature – we somewhat surprisingly found that the relic abundance predicted in the standard approach can differ by up to an order of magnitude from the correct treatment presented in this paper. This is rather remarkable not only in view of the simplicity of this model, but also because the affected region in parameter space happens to coincide with the best-fit region resulting from most recent global scans. We thus expect our results to have a noticeable phenomenological impact, and that our treatment will prove useful also when applied to other examples of relic density calculations in cases where the standard assumption of local thermal equilibrium during freeze-out is not exactly satisfied.

Note added. While preparing this work, we became aware of a dedicated study on resonant DM annihilation Duch and Grzadkowski (2017), which also found that DM can kinetically decouple much earlier than usual in this case.

Acknowledgments.—

We thank Mateusz Duch, Joakim Edsjö, Bohdan Grza̧dkowski, Andreas Hohenegger and Ayuki Kamada for very useful conversations during the preparation of this work. We are also grateful to Tomohiro Abe for pointing out a previously missing factor of 4 in Eqs. (42,43), as well as our detailed discussions around the coupled system of Boltzmann equations (27,II.2). AH is supported by the University of Oslo through the Strategic Dark Matter Initiative (SDI). MG and T. Binder have received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 690575 and No 674896. T. Binder gratefully acknowledges financial support from the German Science Foundation (DFG RTG 1493).

Appendix A Phase-space density evolution of the Scalar Singlet

In Section III.2, we investigated the impact of our improved treatment of the Boltzmann equation on the expected DM relic abundance in the Scalar Singlet model. Here, we supplement this by discussing in some more detail the evolution of the DM phase-space density. The main focus of this discussion, however, will be a more thorough qualitative understanding of the specific features seen in Fig. 2 and Fig. 3, and the underlying interplay of chemical and early kinetic decoupling. Specifically, we can distinguish three mass regimes:

A regime with $53\,{\rm GeV}\lesssim m_{S}\lesssim 60.5$ GeV, which we will refer to as *sub-resonant * because $f_{\chi}$ starts to deviate from its equilibrium value, $f_{\chi,{\rm eq}}=\exp(-E/T)$, at a temperature where the typical DM momenta are too small to hit the resonance, i.e. $\sqrt{s}\lesssim m_{h}-\Gamma_{h}$. As a result, we have $\langle\sigma v\rangle_{\rm(neq)}<\langle\sigma v\rangle_{2,{\rm(neq)}}$333For the sake of better readability, we will suppress the subscript ‘$\rm neq$’ for the remainder of this section. We note that, once chemical decoupling has started, the contribution of thermal averages without this subscript is suppressed by a factor of $Y_{\rm eq}/Y$ in Eqs. (27, II.2).

during the whole freeze-out process in this regime — this is because $p^{2}f_{\chi}(p)$ peaks at a higher value of $p$ than $f_{\chi}(p)$, which brings its bulk distribution closer to (or even on) the cross-section resonance. 2. 2.

A regime with $60.5\,{\rm GeV}\lesssim m_{S}\lesssim 62.5$ GeV that we will refer to as resonant. Here, we have $\langle\sigma v\rangle>\langle\sigma v\rangle_{2}$ around the time when the DM particles start to leave thermal equilibrium, because the larger mass combines with the relevant momenta to $s\sim m_{h}^{2}$. At slightly later times, on the other hand, still relevant in changing the DM abundance, the DM momenta have redshifted so much that we are back to a situation where typically $\sqrt{s}\lesssim m_{h}-\Gamma_{h}$ and hence $\langle\sigma v\rangle<\langle\sigma v\rangle_{2}$. 3. 3.

Finally, there is a super-resonant regime with $62.5\,{\rm GeV}\lesssim m_{S}\lesssim 65$ GeV, where decoupling occurs at such high temperatures that we have $\langle\sigma v\rangle>\langle\sigma v\rangle_{2}$ during the whole time it takes for $Y(x)$ to reach its asymptotic value (determining the relic density).

To help our discussion, let us look at a selection of benchmark points with Scalar Singlet masses $m_{S}=45,57,60.5,62,62.5,63$ GeV and coupling constants $\lambda_{S}(m_{S})$ that result in the correct relic density in the standard approach (dotted line in Fig. 1). In Fig. 5, we show the DM distribution function for these benchmark points that we find with our full numerical approach, for selected values of $x$, and in Fig. 6

the full evolution of $Y(x)$ and $y(x)$ for the different approaches. These figures thus extend the information in Fig. 4 by covering a range of DM masses.

The first thing to note, as exemplified by the benchmark points with $m_{S}=45$ GeV and $m_{S}=63$ GeV, is that for masses sufficiently far away from the resonance we find a phase-space distribution which remains almost exactly Maxwellian in shape. For these points, we therefore find as expected a very good agreement for the evolution of $Y(x)$ and $y(x)$ when comparing the numerical solution and the coupled Boltzmann approach, as well as with $Y$ in the standard Gondolo & Gelmini setup (which assumes $T=T_{\chi}$). We note that this provides an important consistency check for both methods.

An example for a model in the sub-resonant region is the case with $m_{S}=57$ GeV, which we discussed in the main text. Here, the resonant annihilation depletes $f_{\chi}(q)$ for momenta just above the peak of the distribution, leading to a relative decrease with respect to a thermal distribution at these momenta, and hence a decrease in the DM velocity dispersion (aka ‘temperature’). This effect is visible in Fig. 5 starting with a slight suppression at $q\sim 8$ for the curve with $x=20$ (note that the relative enhancement at larger values of $q$ is not relevant for our discussion given that $f_{\chi}$ is already highly suppressed here), and results in the decrease in the evolution of $y$ seen in Fig. 6. The latter can also directly be understood from inspection of Eq. (II.2): in the sub-resonant regime we have $\langle\sigma v\rangle<\langle\sigma v\rangle_{2}$, which drives $y$ to smaller values after decoupling (with a strength proportional to $Y$ – which explains why the scattering term $\propto\gamma(T)$ can increase $y$ again, slightly, once the DM abundance has decreased sufficiently). A second effect of this depletion in $f_{\chi}(q)$ is that $\langle\sigma v\rangle_{\mathrm{neq}}$ decreases, which in turn leads to an earlier chemical decoupling and hence an increased relic density. The difference between the numerical and the coupled Boltzmann approach can in this case thus exclusively be understood as resulting from the slight offset in the $y(x)$ curves during the freeze-out (which in turn results from the fact that the scattering term is not strong enough to maintain an exact Maxwellian shape of $f_{\chi}(q)$ when the velocity dispersion decreases as explained above.)

As we increase the DM mass, we leave the sub-resonant regime and enter the resonant regime, with the transition point marked by the benchmark model with $m_{S}=60.5$ GeV. We note that this transition is also clearly visible in Fig. 2, as a sharp decrease in the temperature at which the DM velocity dispersion deviates from its equilibrium value. The origin of this feature is not an actual delay of kinetic decoupling, but that DM annihilation now starts to deplete $f_{\chi}(q)$ below the peak of the would-be Maxwellian distribution.444 In a similar way, the sharp rise around $m_{S}\sim 54$ GeV in Fig. 2 should not be interpreted as a feature in the momentum exchange rate $\gamma(T)$. Rather, it can be understood as the point where the shape of the $y(x)$ evolution starts to develop from something close to the one in the top left panel in Fig. 6 into something that is much closer to the one in the top center panel (which in turn is driven by the annihilation terms, as explained in the text). As a result, the temperature at which $y$ departs from $y_{\mathrm{eq}}$ increases very quickly as the mass increases beyond this transition point.

This leads to an increase of the velocity dispersion, once equilibrium is left, rather than a decrease as in the sub-resonant regime. This effect is very clearly seen in Figs. 5 and 6, up to DM masses at the higher end of this regime, where the influence of the resonance starts to become less important because we have $\sqrt{s}\lesssim m_{h}+\Gamma_{h}$ only for DM momenta well below the peak of the phase-space distribution.

In the super-resonant regime with $m_{S}\gtrsim m_{h}/2$, finally, we have necessarily $\sqrt{s}\gtrsim m_{h}$. A resonantly enhanced annihilation rate is thus only possible for a very small portion of phase-space, with almost vanishing relative DM momenta. This implies not only that we always have $\langle\sigma v\rangle>\langle\sigma v\rangle_{2}$ in this regime, but also that the effect of the resonance rapidly becomes negligible.

Lastly, it is interesting to note that for $\sqrt{s}\gtrsim m_{h}$ the annihilation rate effectively features a $1/v^{2}$ velocity dependence. This is similar to resonant Sommerfeld-enhanced annihilation, which leads to a suppressed relic density after a prolonged freeze-out phase van den Aarssen et al. (2012). This can clearly be seen in the evolution of $Y(x)$ in Fig. 6, for $m_{S}\sim m_{h}/2$, where the differences between the numerical and the coupled Boltzmann approach are mostly due to the late-time differences in $y(x)$ – which in turn come about because of the rather significant differences in $f_{\chi}(q)$ at large values of $x$ (c.f. Fig. 5).

Appendix B Semi-relativistic kinetic theory

In this Appendix, we discuss how to generalize the highly non-relativistic elastic scattering term in Eq. (5) to incorporate the most important relativistic corrections needed for the numerical implementation of the full Boltzmann equation. Throughout, we refer to this result as ‘semi-relativistic’ scattering.

The starting point is to expand the full collision term $C_{\text{el}}$ in small momentum transfer compared to the typical DM momentum – similar to what is done in order to arrive at Eq. (5), but not only keeping lowest-order terms in $\mathbf{p}^{2}/m_{\chi}^{2}\sim T/m_{\chi}$. From this, we can derive a Fokker-Planck scattering operator in a relativistic form (for details, see Binder et al. (2016)):

[TABLE]

Being a total divergence, this scattering operator manifestly respects number conservation, as it should. Another important property, which one can directly read off from the part inside the brackets, is that it features a stationary point given by the relativistic Maxwell-Boltzmann distribution,

[TABLE]

The non-relativistic limit of Eq. (44) gives the scattering operator (5), but in this limit the stationary point would instead be the non-relativistic version $f_{\chi}^{\text{eq}}\propto\exp[-p^{2}/(2m_{\chi}T)]$ — which would cause a problem in the full BE as this does not correspond to the actual equilibrium distribution fed into the annihilation term of Eq. (37).

In general, the momentum transfer rate $\gamma(T,\mathbf{p})$ in Eq. (44) depends on the DM momentum $\mathbf{p}$. However, the stationary point is independent of $\gamma$, which motivates us to restrict ourselves to the leading order term $\gamma(T)\equiv\gamma(T,\mathbf{0})$, neglecting any momentum dependence, and use the non-relativistic limit in Eq. (44) only to evaluate the momentum transfer rate $\gamma(T)$ as it appears in Eq. (6). To this order, we could thus also replace the leading $E$ in Eq. (44) by $m_{\chi}$; here, we choose to still keep it as it leads to a much more compact analytical form of the equation governing the DM temperature (see below). Explicitly performing the first partial derivative in $C_{\text{el}}$ then leads to the final form of our semi-relativistic Fokker-Planck operator as given by Eq. (8). This operator is our default choice for the numerical implementation of the full Boltzmann equation.

As already pointed out in Section II.3, it is mandatory for the full phase-space calculation to have a scattering operator with a fixpoint that matches the equilibrium distribution of Eq. (45) assumed in the annihilation term. For the coupled integrated Boltzmann system, on the other hand, this issue is fully addressed by using the relativistic temperature definition of Eq. (21) — rather than its non relativistic version typically adopted in the literature in the context of kinetic decoupling — because this automatically leads to the correct fixpoint $T_{\chi}=T$ for both the semi-relativistic Eq. (8) and, to the lowest order, for the non-relativistic version Eq. (5); see the discussion in Section II.2.

Another advantage of our semi-relativistic Fokker-Planck operator is that the differential equation for $T_{\chi}$, often quoted when discussing kinetic decoupling, takes a very simple form even beyond the highly non-relativistic limit. To see this, let us for the moment ignore the impact of annihilations, and take the second moment of the Boltzmann equation with this operator (using the relativistic definition of $T_{\chi}$). This leads to

[TABLE]

which of course is equivalent to Eq. (II.2) in the main text, when neglecting the annihilation terms and implementing the replacement given in Eq. (35). Let us repeat that the r.h.s. of the above equation only takes this particular form with our default choice of the semi-relativistic Fokker-Planck term, whereas the moment appearing on the left hand side is an exact result. This equation is in general not closed in terms of $T_{\chi}$. However, if we make the ansatz of a Maxwellian DM phase-space distribution, c.f. Eq. (34), we get a relation between the different momentum moments,

[TABLE]

such that the differential equation closes in terms of $T_{\chi}$. Indeed, introducing

[TABLE]

it takes a very simple form:

[TABLE]

This generalizes the highly nonrelativistic result Bringmann and Hofmann (2007), for which $w\to 1$ and we hence find the familiar scaling $T_{\chi}\propto T^{2}$ after kinetic decoupling (i.e. when $\gamma\ll H$). In the ultra-relativistic limit, on the other hand, we have $w(T_{\chi})\rightarrow 1/2$ and the likewise familiar scaling of $T_{\chi}\propto T$ for relativistic particles. We note that in the region $x\gtrsim 10$ relevant for early kinetic decoupling, the correction to the non-relativistic limit is already sizeable; e.g. $w(x=10)\approx 0.8$.