DRAKE: Dark matter Relic Abundance beyond Kinetic Equilibrium

We introduce DRAKE, a numerical precision tool for predicting the dark matter relic abundance also in situations where the standard assumption of kinetic equilibrium during the freeze-out process may not be satisfied. DRAKE comes with a set of three dedicated Boltzmann equation solvers that implement, respectively, the traditionally adopted equation for the dark matter number density, fluid-like equations that couple the evolution of number density and velocity dispersion, and a full numerical evolution of the phase-space distribution. We review the general motivation for these approaches and, for illustration, highlight three concrete classes of models where kinetic and chemical decoupling are intertwined in a way that quantitatively impacts the relic density: i) dark matter annihilation via a narrow resonance, ii) Sommerfeld-enhanced annihilation and iii) `forbidden' annihilation to final states that are kinematically inaccessible at threshold. We discuss all these cases in some detail, demonstrating that the commonly adopted, traditional treatment can result in an estimate of the relic density that is wrong by up to an order of magnitude. The public release of DRAKE, along with several examples of how to calculate the relic density in concrete models, is provided at drake.hepforge.org


Introduction
The most often studied scenario to explain the observed cosmological dark matter (DM) abundance [1] considers new elementary particles that have been thermally produced in the early universe. This freeze-out scenario [2] provides an intriguing solution to the DM puzzle not only for classical Weakly Interacting Massive Particle (WIMP) candidates [2][3][4][5][6], with masses at the electroweak scale, but also for much lighter DM particles (which is sometimes referred to as 'WIMP-less miracle' [7]). Various numerical codes are available to provide precision calculations of the DM relic density in these models, including public tools like DarkSUSY [8], MadDM [9], micrOMEGAs [10] and SuperISOrelic [11]. In fact, precision calculations are necessary in order to match the percent-level observational accuracy. In global fits of the full underlying model parameter space, for example, the relic density often provides one of the most relevant observables in terms of contributions to the total likelihood [12][13][14][15][16][17][18][19][20].
All these codes currently rely on the standard approach [21,22] to calculate the relic density, 1 which can be re-cast in the form of a single Boltzmann equation for the DM number density even when including sharp resonances, thresholds or co-annihilations (all of which were initially considered complications [23]). One of the main underlying assumptions for this formulation is that DM (along with potential co-annihilating particles) is kept in kinetic equilibrium with the standard model (SM) heat bath during the entire chemical decoupling process. There exist however various well-motivated scenarios where this assumption is not satisfied, at least not for all relevant parts of the parameter space, and where the standard approach is thus not directly applicable. A possible solution that has received increased attention is to derive coupled equations for the number density and higher moments of the DM phase-space distribution, making use of the Boltzmann hierarchy [24][25][26][27][28][29]. More recently, there have also been attempts to solve the full Boltzmann equation directly at the phase-space level [30][31][32].
Here we introduce a new public code, DRAKE, that is written in Wolfram Language 2 and implements both of these more general approaches for the case of annihilating DM. DRAKE thus complements existing codes to calculate the relic density also in situations where the underlying assumptions of the traditional approach are not satisfied. Additionally, it allows in fact to examine the validity of these assumptions explicitly. As an application, we also study in some detail three concrete physics scenarios where kinetic decoupling can interfere with the freeze-out process: i) s-channel resonances, ii) Sommerfeld enhancement and iii) sub-threshold annihilation. In all these cases, we contrast the results of classical relic density calculations with the more accurate results obtained by DRAKE.
This article is organized as follows. We start in Section 2 by introducing the relevant Boltzmann equations to describe the interaction of DM particles with the thermal bath. In Section 3, the structure of the code is introduced, including the main implemented algorithms that are used to solve these Boltzmann equations. We then discuss in some detail the results of relic density calculations in concrete physics applications in Section 4, before concluding in Section 5. In two Appendices we provide a quick-start guide for using DRAKE (App. A) and discuss examples of elastic scattering operators that can be treated beyond the commonly adopted Fokker-Planck approximation (App. B).

Full Boltzmann equation
We will consider situations where DM interactions with the SM heat bath, through elastic scattering and annihilation processes, are initially strong enough to establish both chemical and kinetic equilibrium. As the universe expands, the DM particles (denoted by χ) fall out of equilibrium and eventually establish the present relic abundance. The evolution of the DM phase-space density f χ (t, p) during this process is governed by what we will refer to as the full Boltzmann equation (fBE): where (2) χχ→f f f χ (E)f χ (Ẽ) describes the effect of two-body annihilations, and the collision term for elastic scattering processes is given by ×(2π) 4 δ (4) (p +k − p − k)|M| 2 χf ↔χf × 1 ∓ g ± (ω) g ± (ω)f χ (Ẽ) − (ω ↔ω, E ↔Ẽ) . 3 In the above expressions, H =ȧ/a is the Hubble parameter, a the scale factor, and we have assumed a Friedman-Robertson-Walker universe, such that f χ only depends on the absolute value of the DM momentum, p = |p|. Furthermore, both collision terms and the squared amplitudes |M| 2 for the respective process are implicitly summed over all heat bath particles f , and final and initial state internal degrees of freedom, respectively. The phase-space distribution of the heat bath particles is given by the usual g ± (ω) = 1/ [exp(ω/T ) ± 1]. Since we assume that DM is non-relativistic around freeze-out, we have neglected Bose enhancement and Pauli blocking factors for f χ (which implies that these factors are also negligible for the heat bath particles in C ann due to energy conservation). For further details about Eq. (1), see Refs. [33,34].

Evaluating the collision terms
In practice, the greatest obstacle in solving Eq. (1) to sufficient precision is often the evaluation of the collision integrals on the right-hand side. For the annihilation term, this is somewhat less critical since -under the generic assumptions of CP invariance and f χ 1 -the phase-space integrals can always be reduced to only one remaining angular integration, and one in energy [21]: where the non-relativistic DM particles in equilibrium follow a Maxwell-Boltzmann distribution, f χ,eq (E) = e −E/T . Further, σ is the annihilation cross-section for the processχχ →f f , and v M ≡ [s(s − 4m 2 χ )] 1/2 /(2EẼ) is the Møller velocity; in the rest frame of either of the DM particles, it equals the velocity of the other particle, v M = v lab ≡ [s(s − 4m 2 χ )] 1/2 /(s − 2m 2 χ ). For later reference, let us here also introduce the angular (Ω) averaged quantity σv θ ≡ 1 2 d cos θ v M σχ χ→f f that is used in our numerical code and depends only on the magnitude of the three-momenta, p andp.
The elastic scattering term also simplifies for highly non-relativistic DM. For m f m χ , specifically, the typical momentum transfer per collision is then much smaller than the average DM momentum in equilibrium. Expanding C el up to second order in the momentum transfer results in a simple differential operator of Fokker-Planck type [35][36][37], which can be used to describe kinetic decoupling taking place much later than chemical decoupling (see, e.g., Ref. [38]). To improve the description of early kinetic decoupling, we keep in DRAKE some of the leading relativistic corrections [30,39], resulting in the Fokker-Planck operator (FP): It is important to note that C FP [f χ ] = 0 for a relativistic Maxwellian shape f χ ∝ e −E/T , which is consistent with the stationary solution of the relativistic annihilation term in Eq. (4) (see appendix in Ref. [30] for more details). Furthermore, it is worth mentioning that the above operator can be written as a total momentum divergence and hence manifestly conserves the DM particle number. The momentum transfer rate γ(T ) introduced above is the same as in the highly non-relativistic version of Eq. (5), and given by (see also Refs. [40,41]) where k 2 cm ≡ m 2 χ k 2 /(m 2 χ + 2ωm χ + m 2 f ) and the scattering amplitude entering the differential cross-section, , is evaluated at s m 2 χ + 2ωm χ + m 2 f . We will encounter one concrete example in this work, in Section 4.3, where the mass of the scattering partner is comparable to the DM mass, such that the momentum transfer can no longer be assumed to be small. In this case, we have to resort to the full scattering term in Eq. (3), which we re-write in a form that is more suitable for numerical integration, and which formally resembles the annihilation case, Eq. (4): 3 Here, the quantity W is defined as and has the same physical dimension as a cross-section. For our numerical codes it is convenient to factorise out its angular average W Ω ≡ 1 4π dΩ W , depending only on T , p andp. This quantity is in general challenging to compute due to the multidimensional integrals. In Appendix B we summarize, however, some concrete cases (including the one relevant for the discussion in Section 4.3) where the amplitude M takes a form allowing this general expression for C el [f ] to be analytically reduced to a one-dimensional integral.

Fluid equations
An alternative to the numerically challenging task of directly computing the evolution of the phase-space density f χ consists in restricting the discussion to its first moments. Instead of the full Boltzmann Eq. (1) one thus only considers the corresponding momentum moments of this equation, thereby implementing a 'hydrodynamical' approach to the problem that essentially results in a set of fluid equations. Just as in the case of hydrodynamics, however, additional assumptions are needed to close the Boltzmann hierarchy at any given level.
The simplest case, and in fact the traditional way to handle relic density calculations, is to only consider the lowest moment of f χ , namely the DM number density 3 . The additional requirement to close the resulting equation for n is typically met by assuming that kinetic equilibrium is maintained until the chemical decoupling process is completed. For nonrelativistic DM, this implies that the DM phase-space distribution is of the form f χ = exp[(µ − E)/T ], which fixes f χ up to a function of T (here parameterized as the DM chemical potential µ, in this case related to the number density as µ/T = ln[n/n eq ]).
Inserting this form into Eq. (1), and integrating over the DM momenta p, then results in what we will refer to as the traditional number density equation (nBE): where n eq refers to the DM density in chemical equilibrium, i.e. for µ = 0. The thermal average ... T of the velocity-weighted annihilation cross-section can be stated in terms of a single integral over the center-ofmass energy, as explicitly given in Eq. (3.8) of Ref. [21]. In situations where kinetic decoupling interferes with the chemical decoupling process, the main assumption leading to Eq. (9) is clearly too restrictive. This implies that one needs to move (at least) one level up in the Boltzmann hierarchy to (approximately) describe such situations, and thus to leave the second moment of f χ as a dynamical degree of freedom. A convenient parameterization for this is the DM velocity dispersion or 'temperature', and one way of closing the Boltzmann hierarchy at this level is to assume, in analogy to the case discussed above, f χ = exp[(µ − E)/T χ ] (with µ/T χ = ln[n/n eq (T χ )]; see Ref. [30] for a more detailed discussion). Using Eq. (5), and assuming entropy conservation, this results in what we will refer to as the coupled Boltzmann equations (cBE): parameterizes the deviation from the highly non-relativistic case (w = 1). In order to arrive at this compact form, we have followed the standard convention of introducing dimensionless variables Y (x) ≡ n/s, Y eq (x) ≡ n eq (T )/s, where s is the entropy density. 4 The thermal average σv 2,T is a variant of the commonly used thermal average σv T , and is explicitly stated in Ref. [30]. We also introduced which in the Fokker-Planck approximation of Eq.
, with g s eff being the entropy degrees of freedom of the background plasma. The first of these equations, Eq. (10), is the analogue of the traditional number density equation; in the limit T χ → T (kinetic equilibrium) it reduces as expected exactly to Eq. (9) expressed in dimensionless variables.
Let us close this section by briefly mentioning the potential effect of DM self-interactions [42,43] on our discussion. In principle, this would be described by an additional collision term C self on the right-hand side of Eq. (1), but bringing such a term into a form that is numerically tractable is challenging. In light of this situation it is interesting to note that the two beyond-nBE approaches implemented in DRAKE can be regarded as an effective way of bracketing the impact of such model-dependent DM self-interactions: our implementation of fBE provides the correct description of how the DM phase-space density evolves if the effect of DM self-interactions is negligible, while the Eqs. (10,11) that cBE is based on become exact under the assumption that DM self-interactions are maximally efficient and hence force the DM distribution to be of the form

Code design
DRAKE is a package of routines written in Wolfram Language that allows to numerically compute the DM relic abundance, Ω χ h 2 , in the nBE, cBE and fBE approaches introduced above. In Section 3.1, we provide an overview of the code's modular structure and how to use it in practice to obtain Ω χ h 2 and other related quantities for a given DM particle model. The internal algorithm implemented for the time (x) integration of the three Boltzmann equations is presented in Section 3.2, and in Section 3.3 we summarize various validity checks of the code that we have performed.

Overview
The DRAKE package provides routines to solve (time integrate) the individual Boltzmann equations, as well as a set of preparatory routines that compute quantities required by the solvers, e.g., averages of the annihilation cross-section. This design structure allows to perform relic abundance computations in a flexible manner, to reuse certain output quantities from preparing routines in multiple Boltzmann solvers, and reduce the required model input to fairly simple expressions.
To make the code design concrete, let us start with the three Boltzmann solver routines nBE, cBE and fBE listed in Table 1. As explained in Section 3.2, these are adaptive and implicit ordinary differential equation solvers, where the partial differential Eq. (1) for the fBE solver has been mapped to a set of coupled ordinary differential equations by the so-called method of lines (see, e.g., Ref. [44]). Each solver calculates the relic abundance Ω χ h 2 and the full Y (x) evolution. cBE and fBE also give the y(x) evolution, and the latter additionally provides information on the full phase-space density f χ (x, p).
The annihilation cross-section averages (third column) needed in the Boltzmann equations are computed by the preparatory routines PrepANN, PrepANN2, and PrepANNtheta, as listed in Table 2. The only input required from the user is the DM annihilation cross-section in the form of σv lab as a function of the Mandelstam variable s. The routines then adopt to adjustable accuracy settings and the DM model at hand, in order to generate accurate and densely tabulated outputs for the resulting interpolation functions.
For the cBE and fBE routines also the scattering operator C el in Eq.    The preparatory routines, computing averages of the velocity-weighted annihilation cross-section σv lab . Output of these routine calls is stored in the session memory, such that, e.g., the quantity σv can be used for both nBE and cBE routine calls. These routines are located in rates.wl (in the src/ directory).
amplitude |M| 2 χf ↔χf -should be provided by the user. The preparatory routine PrepSCATT then adaptively tabulates γ(x) to provide an interpolating function for fast evaluation. However, we stress that the cBE and fBE routines are not limited to the use of the Fokker-Planck approximation for the scattering collision term C el . As explained in Appendix B, the full C el can be used if the quantity W Ω from Eq. (8) is provided (which in practice amounts to providing C el 2 andĈ el as defined in B.2 for cBE and fBE, respectively).
The first code release also comes with five preimplemented DM model setups. These are the Scalar Singlet DM model [45][46][47] 5 , our three example scenarios discussed in Section 4, and a WIMP-like toy model. Each of these comes with three files: a model file (containing σv lab (s) and either γ(x) or the scattering operators C el 2 andĈ el ), a parameters file (with numerical values of the DM mass, couplings, internal d.o.f., etc.) and a settings file (with various code options).
The DM relic abundance in these, or any other user-defined, setups can be conveniently calculated, e.g. within a Mathematica notebook, by sequentially calling the required DRAKE routines. A practical alternative is to use the customizable template script main.wls provided in the main directory. As illustrated in Fig. 1   . The command GetModel then loads the content from these <files> before the script continues to execute the required preparatory routines. The preparatory routines PrepANNtheta and PrepSCATT will not be called if the options Analyticsvtheta or FullCell are set to True, respectively; in that case the quantities σv θ or, respectively, C el and C el 2 , need to be provided by (any of) the user defined <files> (as indicated by the green texts along the dashed lines).
The main physics output of each preparatory routine is stated in the respective red box right below, and the subsequent calls of the Boltzmann equation solvers nBE, cBE and fBE depend on these physical quantities as indicated by their incoming arrows. The output from each Boltzmann solver, finally, is explicitly stated in the red boxes in the bottom row.
calculation procedure consists of initializing the DRAKE package, then loading the DM model setup files, performing the required preparatory calculations and finally executing the Boltzmann equation solvers.
A dedicated quick start guide as well as a description of more features and details of DRAKE are provided in Appendix A.

Numerical implementation
All three Boltzmann equations (nBE, cBE and fBE) are numerically time (x) integrated by essentially the same implicit adaptive algorithm. Once the phasespace density is discretized in momentum space as in particular, all of them can be written as a coupled system of ordinary first order differential equations of the form where for cBE as in Eqs. (10) and (11), and V(x) = {f 0 (x), ...f N (x)} for the momentum-discretized phasespace density in the fBE approach. We follow standard practice for the numerical integration of these equations, see, e.g., Ref. [48]. Concretely, we choose the Adams-Moulton time discretization methods of order one and two, which are also known as the implicit Euler and implicit trapezoidal method where The adaptive step size h is controlled through a local relative error of the two methods,  (14) and (15) can respectively be solved analytically for V T i and V E i , allowing for efficient time integration. In this one-dimensional case, thus, the implemented algorithm in our nBE routine is equal to the one used in DarkSUSY.
In the case of cBE and fBE, the solution of Eqs. (14,15) needs to be obtained numerically. A standard root finding algorithm is the Newton iteration method. Defining the iteration depth n and introducing , the Newton iteration for the Euler method in Eq. (14) can be written in the form of a linear algebraic problem for ∆ , and for the trapezoidal method in Eq. (15) as where ∂ V F denotes the (N + 1 × N + 1) Jacobian matrix. The improved value for each iteration can be obtained as V . At each step in the x integration, we choose the predictor as V (0) i = V i−1 and iterate each method until the maximum relative error, errNewton is one order of magnitude smaller than the upper bound on err i . We also require errNewton (n+1) i to be smaller or equal to another accuracy control parameter, errMaxNewton i = 0.4, at every iteration step; otherwise, the stepsize h is lowered in order to ensure that the Newton iteration converges. Default values for the accuracy control parameters introduced above (err i , errNewton (n+1) i and errMaxNewton i ) can be adopted by changes in the settings file, see Appendix A.5.
For the cBE and fBE routines we have implemented further, automatized performance optimizations. In the two-dimensional cBE case, in particular, we analytically solve in Eqs. (14,15) the corresponding number density equation for Y i , Eq. (10), as a function of the DM temperature y i , and plug that solution into the DM temperature equation for y i , Eq. (11). Effectively, this reduces the cBE system to an one-dimensional ordinary differential equation, which is how it is implemented in the cBE routine.
For fBE, the computation of the Jacobian ∂ V F is implemented in vectorized form. Moreover, all parts in F that depend only linearly on V, e.g., the scattering operator in Eq. (5), are pre-computed for every time step (i) and separated from the ones which need to be updated in every Newton step (n).
Furthermore, we implemented two momentum coordinate systems, A and B, for different temperature regimes, such that the phase-space density long before and long after kinetic decoupling remains stationary in the respective coordinates (for Y = const.). Concretely, the phase-space density is disretized on a uniform grid in q A ≡ p/ m χ T for Abs[1 − y eq /y] <qATOqB and q B ≡ (g s eff ) −1/3 p/T otherwise. By default, this setting variable is set to qATOqB= 0.1.
For the momentum derivatives of the phase-space density in Eq. (5) we implemented finite differentiation methods of higher order accuracy using several neighboring points. For central differentiation and fourth order accuracy implemented as a default, the boundary values of the phase-space density used are reflection symmetry at the origin f 0 , introducing the ghost points to obey , and vanishing phase-space density at numerical infinity, i.e., In practice we rely on Mathematica's LinearSolve command for solving Eqs. (16) and (17) in the fBE approach. This part of the code, as well as other timeconsuming matrix manipulations like collision term updates, is implemented in a C-compiled environment (if no C compiler is available it is targeted to the Wolfram Virtual Machine), allowing for run-time performances comparable to pure C++ or Fortran code based on the analogue dgesv command from the LAPACK package [49]. 8

Validity tests
The evolution of the number density Y (x) based on the traditional nBE approach were cross-checked against results obtained with DarkSUSY, for a variety of models ranging from rather generic WIMP realizations to situations where DM annihilates through a narrow resonance. Similarly, the temperature evolution y(x) computed with the kinetic decoupling (only) routines in DarkSUSY was compared to that obtained in DRAKE in the cBE and fBE approaches after switching off annihilation in our code through the settings command KDonly=True. For all tested models the results were found to be in agreement at well below the percent level.
We further checked explicitly that the number density evolution in the cBE and fBE approaches correctly coincides with the nBE solution in situations where the effect of kinetic decoupling is irrelevant; e.g. when the DM particles remain in full kinetic equilibrium during the freeze-out process or when the velocity dependence of the annihilation cross-section is insignificant.
The non-trivial interference of chemical and kinetic decoupling, as taken into account in our two beyond-nBE approaches, was meticulously checked against independently implemented codes. In particular, the DRAKE cBE results in the current implementation of the Scalar Singlet model (c.f. footnote 4) were compared to results obtained from FORTRAN routines relying on the sfode [50] solver, finding agreement at the sub-percent level. 6 Furthermore, the fBE results were compared to the phase-space solver originally developed in MatLab based on the ODE15s solver and used in Ref. [30]. 7 For resonance annihilation considered in Section 4.1 for example, very good agreement was found. Note that the cBE and fBE routines in DRAKE do not rely on pre built-in differential equation solvers (see Section 3.2), so all these tests provide an independent validity check.
A final important consistency check of the numerical fBE implementation concerns the question whether the elastic scattering operator in Eq. (5) conserves particle number. While this is manifestly the case in the continuous limit (and hence also in cBE), the discretized version with a finite number and range in momentum can lead to a spurious change of the comoving number density if γ/H ≫ 1. Therefore some care must be taken in the initial high-temperature regime. One way to reduce these purely numerical artifacts is to increase the number of phase-space density elements, which however affects the run-time of the fBE routine quadratically. We have therefore introduced an upper limit on γ/H above which the system is assumed to be in kinetic equilibrium and DRAKE adopts the nBE approach. By default, this upper limit is set to gamcap= 10 5 . This value and the number of phase-space elements, qN, can be adopted in the settings file. 8 For specific problems, optimal accuracy settings for fBE can differ from the default ones. We therefore provide for all examples in Section 4, as well as for the Scalar Singlet Model, suggested settings files. On a standard laptop computer the time to compute the relic density for one parameter point in the fBE approach can range from less than a second to tens of minutes (e.g. in the presence of very narrow resonances). Furthermore, DRAKE contains example benchmark files with settings adjusted to specific cases, illustrating the usage for more challenging models. For further details concerning those test files see Appendix A.3, as well as Appendix A.5 for accuracy control parameters.

Physics scenarios
A velocity-dependent annihilation cross-section is a necessary condition for deviations from the standard computation of the relic abundance, due to the interplay of chemical and kinetic decoupling. In this section we consider three different types of physically well-motivated scenarios where such a velocity-dependence can appear, and where kinetic equilibrium is not necessarily maintained during chemical decoupling. We discuss these scenarios in some detail, both to highlight the nontrivial physics of the freeze-out process in at least parts of the models' parameter space, and to illustrate how DRAKE can be used to study the freeze-out process beyond the simplifying assumptions usually adopted in the literature.

Resonant annihilation
As our first example, we consider cases where the annihilation cross-section has a strong velocity dependence induced by an s-channel resonance. The effect of early kinetic decoupling for such a setup has been studied in detail in Ref. [30] for the Scalar Singlet model [45][46][47], where the almost on-shell particle in the s-channel is the SM Higgs. For Scalar Singlet masses slightly smaller than half of the Higgs mass it was found that a larger coupling to the Higgs is required to obtain the right 9 relic abundance, which interestingly implies that future measurements of the Higgs-to-invisible decay width can probe more of the parameter region than expected from the standard relic abundance computation, see Refs. [25,29,51] for experimental probes and also for other Higgs resonance scenarios with similar effects.
To offer a complementary perspective, we consider here instead a generic vector mediator A µ inducing an s-channel resonance (where the parameter region close to the resonance is particularly interesting from a modelbuilding perspective, see, e.g., Refs. [52][53][54][55]). A minimal version of such a scenario is described by the interaction Lagrangian allowing a Dirac fermion DM particle χ to resonantly annihilate into heat bath fermions f . The model can thus be described by a set of 5 parameters, (m χ , r,γ, δ, ρ), where m χ is the DM mass, r ≡ m f /m χ defines the mass ratio of heat bath fermions to DM,γ ≡ Γ A /m A is a dimensionless measure of the total decay width of A µ , ρ ≡ √ g χ g f is the combination of coupling constants relevant for our discussion and δ ≡ (2m χ /m A ) 2 − 1 measures the deviation from the exact resonance position (at δ = 0).
With these definitions, the annihilation cross-section for this process, χχ → A → ff , can be written as with α(s) = 4(2s + 1)(2s + r 2 ) ands ≡ s/(4m 2 χ ), where √ s is the center-of-mass energy, and the resonance is encoded in the Breit-Wigner propagator 9 The scattering amplitude for χf ↔ χf in the nonrelativistic limit,s → (1 + r 2 + 2ω/m χ )/4, on the other hand, is given by where Following our standard convention, Eq. (21) is summed over both in-and out-going spin states, as well as particle and anti-particle states of the scattering partners f . 9 For gχ g f and |δ| 1, Eq. (20) may require a modification due to an energy-dependent width [25]. We will concentrate here on the parameter region |δ| < 1 where the annihilation cross-section is resonantly enhanced. This implies that ρ must be suppressed in order to obtain the correct relic abundance, leading to a corresponding suppression of the scattering amplitude. On top of that, the momentum transfer rate γ(x) can be further Boltzmann suppressed for sufficiently massive scattering partners. In combination, these effects can lead to very early kinetic decoupling, interfering with the chemical decoupling process. In Fig. 2 we quantify the impact on the relic abundance in such a situation, choosing for definiteness a fixed DM mass of m χ = 1 TeV and a resonance width ofγ = 10 −5 . Here, we fix the coupling ρ in every parameter point such that the conventional approach (nBE) matches the value of the observed DM abundance, i.e., 2(Ω χ h 2 ) nBE = Ω DM h 2 = 0.120. Clearly, the cBE and fBE approaches demonstrate significant corrections to the standard treatment. As expected, the larger the mass ratio r (from dotted to dashed to solid lines), the earlier DM decouples and the stronger the effect becomes. 10 Even for very light heat bath particles, however, the correction to the standard result remains sizeable (for r < 0.1, the resulting curves do not significantly differ from the r = 0.1 case displayed in Fig. 2). The characteristic two-bump feature and the narrow dip around the exact resonance position at δ = 0 are a consequence of different DM cooling and heating effects during the chemical evolution that affect the DM phase-space distribution f χ (x, p). The phenomenology of coupled chemical and kinetic evolutions here is very similar to the Scalar Singlet DM example, and for a detailed discussion of the origin of these features we refer to to Appendix A in Ref. [30].
In Fig. 3, we complement this discussion by showing the impact on the relic abundance when instead varying the resonance widthγ and keeping the mass ratio r fixed. This leads to a structure that is straightforward to relate to what is visible in the previous figure; for example, the two distinctive peak regions in the bottom left corner correspond to the two peaks in Fig. 2 (note however that here we consider a much lighter DM particle, m χ = 1 GeV, than in Fig. 2). For most of the parameter space, a smaller width generally leads to a larger effect, i.e. larger deviations from the standard computation. It is interesting to note that even for widths as large as that of the SM Z-boson the refined prediction of the The coupling strength is set to ρ = 1.13 × 10 −2 in order to satisfy 2(Ωχh 2 ) fBE = Ω DM h 2 . The initial equilibrium phasespace distribution is strongly distorted during chemical and kinetic decoupling, and finally remains in a highly non-thermal shape.
DM abundance can deviate at a level well exceeding the typically quoted observational uncertainty of ∼ 1% in Ω DM h 2 -and hence at a level that would, e.g., affect global fits in a noticeable way.
It is worth noting that for the simple model considered here, not every pair of values (δ,γ) shown in Figs. 2 and 3 may be a consistent choice. Indeed, the minimal contribution to the width, from the interaction Lagrangian in Eq. (18), is given bỹ In Fig. 3 we indicate (with gray shaded regions) the values of (δ,γ) that cannot be satisfied by Eq. (23) when ρ = √ g χ g f is fixed by the relic density condition (either because Eq. (23) would imply a larger value ofγ than required, or because at least one of the couplings would no longer satisfy g χ , g f < √ 4π, thus indicating a breakdown of perturbativity). Let us finally remark that the more narrow the resonance, the more momentum-selective is the annihilation process. This can lead to shapes of the distribution function that strongly differ from thermal ones. As a concrete example, we demonstrate in Fig. 4 that the fBE approach implemented in DRAKE can accurately resolve such non-trivial phase-space evolutions even for extremely narrow resonances (in this example a factor of a few below the Higgs width).

Sommerfeld-enhanced annihilation
As our second example, we consider a case where DM annihilation is Sommerfeld-enhanced due to the presence of a light mediator. Physically, such a light mediator induces a long-range Yukawa potential between the DM particles that modifies their wave function, leading to a non-perturbative enhancement of the tree-level annihilation rate [56][57][58][59][60][61]. Because the Sommerfeld effect is strongly velocity-dependent, it provides a prime example of interest to study the interplay of chemical and kinetic decoupling [24,[62][63][64]. In fact, this is the context in which coupled Boltzmann equations akin to our Eqs. (10,11) have first been proposed [24].
For simplicity, we consider the same model Lagrangian as in the previous section, Eq. (18), but now in a very different parameter region. Concretely, we will assume m A α χ m χ , with α χ ≡ g 2 χ /(4π), such that the vector mediator induces a long-range force. Furthermore, we take the heat bath particles f to be (approximately) massless and only milli-charged under the broken U (1) , in the sense that g f g χ . Consequently, χχ → AA is the leading annihilation process. The corresponding tree-level cross-section, expanded up to leading order in v 1 and m A /m χ 1, is given by (σv lab ) 0 πα 2 χ /m 2 χ . We explicitly checked that this approximation holds for the range of parameters that we will consider here. The full s-wave annihilation cross-section including the Sommerfeld effect is then obtained as where we use the analytical expressions for the Sommerfeld enhancement factor S(v lab ) that were obtained, e.g., in Ref. [65,66] after approximating the Yukawa potential by a Hulthén potential: where v ≡ v lab /(2α χ ) and A ≡ m A /(α χ m χ ). We concentrate on the parameter region where g f is (just) large enough for the mediator to be in equilibrium during freeze-out, established by the decay and inverse decay processes A ↔ ff . Requiring the decay rate to satisfy Γ A→ff /H 1 for T m χ , one finds that this is realized for g f 3×10 −6 (0.1/r A )(m χ /TeV) 1/2 , 11 where r A ≡ m A /m χ . For simplicity we will exclusively focus on coupling values that are not significantly larger than this lower bound, thus implementing the earliest possibility when DM can kinetically decouple in this model. In this case, Coulomb-like scattering processes f χ ↔ f χ can be neglected -as we checked explicitly -and DM is kept in kinetic equilibrium by the Compton-like scattering with the vector particle, χA ↔ χA. The corresponding scattering amplitude is for non-relativistic DM given by where, to leading order in the DM mass, . (27) In this parameter region, the phenomenology of the model thus only depends on the coupling g χ (unlike in the previous section, where the relevant coupling combination was ρ = √ g χ g f ). For relevant coupling values and the scattering amplitude as in Eq. (26), kinetic decoupling occurs once the vector boson A enters the non-relativistic regime, causing a Boltzmann suppression of the momentum transfer rate γ(x).
In Fig. 5 we show the resulting relic density in this setup, for a range of mediator masses m A and fixed DM mass (m χ = 2 TeV) and coupling (α χ = 0.07). As is clearly visible in this figure, the Sommerfeld effect alone affects the relic density by a factor of 2-3 compared to the tree-level expectation. For selected values of mediator masses, the difference between the standard approach and both our beyond-nBE approaches can be even larger. For these parameter combinations, bound states with zero binding energy can form; parametrically, these exist precisely at threshold for m A /(α χ m χ ) = 6/(n 2 π 2 ) in the case of the Hulthén potential (with n ∈ Z + ). Close to these parametric resonances the Sommerfeld effect scales as S(v lab ) ∝ 1/v 2 lab for v lab m A /m χ , leading to a second period of annihilation after kinetic decoupling [24,[62][63][64]. We illustrate this in Fig. 6 by plotting the evolution of the DM abundance Y (x) for three selected mediator masses close to such a parametric resonance (left panel). Due to the inverse velocity dependence of S(v lab ), furthermore, DM particles prefer 11 Since m A mχ ∼ T we need to take into account relativistic corrections to the decay rate in the form of a time dilation factor m A /E, see, e.g., Ref. [67]. We estimate the decay rate for m A T as Γ A→ff is the decay rate in the A rest frame and m A /(2T ) is the average time dilation in the plasma frame. In this estimate, Pauli blocking of final states was neglected and we adopted a classical Maxwell-Boltzmann statistics for the vector boson A. to annihilate at smaller momenta, leading to a selfheating after kinetic decoupling and hence an increase in y(x) (right panel) -see also the discussion in Ref. [24].
Inspecting Fig. 5, we find significant deviations from the standard computation also for parameter regions further away from the exact resonance condition. Here, the Sommerfeld effect scales only as S(v lab ) ∝ 1/v lab , which correspondingly leads to a significantly weaker re-annihilation effect. Let us point out that even for situations with a relatively early kinetic decoupling as studied in this section (though still much later than what we studied in Section 4.1), the difference between cBE and fBE is as expected rather small, independently of the vicinity to a parametric resonance. This holds both for the evolution of the number density and that of the DM temperature, as also demonstrated earlier in Ref. [68].
In general, the effect on the DM abundance due to kinetic decoupling is of course larger for earlier kinetic decoupling (assuming all other parameters to be fixed). However, we remark that even for later kinetic decoupling than in the minimal case discussed here -e.g. in the case where Coulomb-like scattering dominates, for much larger coupling values g f -the correct prediction of Ω χ h 2 can be sensitive to the interplay of kinetic and chemical freeze-out processes, making it necessary to go beyond the traditional nBE approach.  Fig. 5. Line colors correspond to nBE (green), cBE (blue) and fBE (red). Line styles refer, as indicated, to different values of m A , chosen such that n 2 = 6αχmχ/(π 2 m A ) is close to integer (parametric resonance condition). The second annihilation era characteristic for these models is clearly visible as a drop of Y in the left panel; the corresponding increase in the DM temperature visible in the right panel is caused by a self-heating due to efficient annihilation out of thermal equilibrium.

Sub-threshold annihilation
As our final example we consider a situation where the annihilation process that sets the relic density is kinematically not accessible in the limit of vanishing velocities, sometimes dubbed 'forbidden' DM [69].
For concreteness, we adopt a simple scenario with an interaction Lagrangian given by where the scalar φ 1 takes the role of the DM particle. It directly interacts only with the scalar φ 2 , which we assume to be close in mass with φ 1 , i.e., r ≡ m 2 /m 1 ∼ 1, and to be in thermal contact with the heat bath fermions f . To lowest order, the total DM annihilation crosssection is determined by the process φ 1 φ 1 → φ 2 φ 2 , and given by For r > 1, this process is only open to particles in the high-energy tail of the DM distribution, leading to an exponential suppression of σv and a corresponding increase of Ω χ h 2 . Similar to the Sommerfeld enhancement example discussed in the previous section, we restrict our discussion for simplicity to couplings y f (just) large enough to bring φ 2 into equilibrium before the onset of the freeze-out 13 process, through (inverse) decays φ 2 ↔f f . Demanding for concreteness Γ φ2↔f f /H 1 for T m 1 , this requires y f 10 −7 r − 1 2 (m 1 /TeV) 1 2 . At the same time we require for consistency that y f is small enough such that i) close to the threshold at r ∼ 1, 3-body processes of the form φ 1 φ 1 → φ 2 φ * 2 → φ 2 ff can be neglected and that ii) the loop-induced scattering with f is subdominant compared to the tree-level scattering with the Boltzmann-suppressed φ 2 particles, via φ 1 φ 2 ↔ φ 1 φ 2 . 12 In this parameter region, the momentum transfer rate, Eq. (5), is then given by We note that this expression explicitly features the already mentioned exponential suppression due to scattering with non-relativistic φ 2 particles. Hence, we expect kinetic decoupling to happen very early in this model, around the time of chemical decoupling. As already stressed in Section 2.2, however, the approximate scattering collision term in Eq. (5), including the momentum transfer rate, relies on the assumption of small momentum transfer compared to the average DM momentum -which is typically not satisfied for scattering partners that are close in mass to the DM particle. For the simple case of a constant scattering amplitude, consistent with our interaction Lagrangian in Eq. (28), we demonstrate in Appendix B that it is possible to instead implement the scattering collision term C el without this assumption. Concretely, we perform analytically all integrals for W in Eq. (8) , c.f. Eq. (B.2), and use this expression to compute C el in Eq. (7). This allows us to contrast our results to those based on the approximate scattering term C FP (relying on small momentum transfer).
In Fig. 7 we plot the relic density that results from the fBE and cBE approaches, with different implementations of the scattering terms, relative to that from the nBE approach. Here, for each value of r, the coupling λ is fixed by the requirement that the standard nBE prediction matches the observed abundance. As is clearly 12 We explicitly checked that these requirements actually allow values of y f up to a few orders of magnitude above the lower bound due to the thermalization condition, with condition i) being the more stringent requirement. Concretely, we followed Ref. [70]  visible, early kinetic decoupling can have a significant impact on the relic abundance also for this threshold example, at least for r 1. While it is remarkable that both scattering term prescriptions show qualitatively the same result, however, in this model DM turns out to be kept slightly more efficiently in kinetic equilibrium with the more correct C el approach; see also Fig. 10 in Appendix B. This explains why the C el prescription leads to a relic density slightly closer to that of the standard nBE approach than what we find with the Fokker-Planck approximation.
In Fig. 8, we show the abundance and temperature evolution for selected values of r > 1. Qualitatively, DM needs higher momenta to overcome the annihilation threshold, leading to a self-cooling phase as soon as it is no longer (fully) kinetically coupled to the particles φ 2 that remain in equilibrium with the heat bath. Due to this cooling, annihilation becomes less efficient earlier, resulting in a higher DM abundance than in the standard computation.
Let us close this section by briefly mentioning again that the kinematically suppressed annihilation crosssection implies that the coupling λ must grow exponentially with r in order to maintain the correct DM abundance. For example, in the nBE approach with m 1 = 100 GeV very large couplings λ 4π are reached at r 1.2; for the cBE and fBE approaches this happens at even smaller values of r, closer to the peak visible in Fig. 7.

Summary
Thermal freeze-out is widely considered one of the most intriguing mechanisms for the production of DM. The underlying assumption of by far most relic density calculations in the literature is that kinetic equilibrium is maintained during the freeze-out. Here we have presented a new public tool, DRAKE, to explicitly gauge the impact of this assumption in cases where it might not be valid. To do so, the code offers various alternative schemes to calculate the relic density that take into account the intriguing effects of kinetic decoupling during the chemical freeze-out process, including a full calculation at the phase-space level.
In fact, it has repeatedly been pointed out that chemical and kinetic decoupling can be intertwined in a way that significantly affects the result of relic density calculations [24][25][26][27][28][29][30][31][32][71][72][73][74][75]. To provide context, we have therefore devoted a large part of this article (Section 4) to a comprehensive and updated study of three physically well-motivated scenarios of annihilating DM where this is the case, illustrating the need to move beyond the standard treatment and demonstrating that this is possible without compromising on accuracy. This is clearly important for global fits that include parameter regions in their scans where the interplay between kinetic and chemical equilibrium cannot be neglected, but can turn out to be relevant also in various model-building considerations.
Let us finally stress that, while the main focus of DRAKE is the relic density, its output is by no means restricted to a single number. Rather, it allows to compute the full time evolution of the DM phase-space density, or its lowest moments, which may be connected to further late-time observables. For example, a non-standard velocity distribution would affect how free streaming impacts the matter power spectrum of density perturbations [76,77], which is one of the main motivations for why linear perturbation solvers like CLASS [78] explicitly support externally tabulated (non-standard) phasespace densities. An exploration of such effects in the context of kinetic decoupling interfering with chemical decoupling would be warranted, but beyond the scope of this work. and documentation of the main functions and variables are provided with the usual syntax, e.g., calling ?cBE returns information about the cBE routine. All input quantities can be loaded with the command GetModel["<model>", "<parameters>", "<settings>"]; The three arguments are names of files (without their extension .wl), located in the sub-directory ./models/< PrepANNtheta σv θ tsvtheta, isvtheta When working with all approaches simultaneously, it is sufficient to call each preparatory routine once. Optional procedures, where e.g. PrepANNtheta and PrepSCATT are omitted, are explained in Section A.5. The output of the preparatory routines and Boltzmann solvers is summarized in Table 3. All listed output variables are stored in the kernel session memory and can be directly accessed and saved to a file <outfile> with, e.g., DRAKE's command save ["<outfile>"]. Values of any model or setting parameters can be changed actively during a notebook session. This makes scans over, e.g., the DM mass simple to perform. Alternatively, one can create different <parameters> and <settings> files and repeatedly use GetModel.

A.2: Template script
The release also contains a template script main.wls, which serves as a streamlined and customizable way to use DRAKE directly from a terminal. For any of the five pre-implemented models (see Section A.1 for file names), in particular, the template script can be executed as Code settings are stored in <settings> files for all preimplemented models. These consist of accuracy control parameters and Boolean variables for optional code usage. A possible way of directly checking the effect of these parameters is to plot the output variables of the relevant routines, as listed in Table 3. All parameters introduced in this section are also briefly described in the <settings> files.
Options The options for preparatory and Boltzmann solver routines are determined by a set of global Boolean variables. In particular, the Boolean variable RelThermalAv determines how σv lab is averaged in the PrepANN and PrepANN2 routines. If set to true, σv lab must be provided as a function of the Mandelstam s variable and averages are computed fully relativistically as in Ref. [21,30]. Setting RelThermalAv to false, σv lab must be provided as a function of v lab , and the calls PrepANN and PrepANN2 compute averages in the highly non-relativistic limit, as in Ref. [24]. The Sommerfeld enhanced annihilation scenario, as described in Section 4.2 and implemented in SE.wl, is one example where σv lab is a sufficiently smooth function for accurate non-relativistic thermal averages. Especially for PrepANN2 this is of great advantage since the amount of numerical integrals to be performed is reduced to one [24]. The Boltzmann solvers cBE and fBE can take optional routines as input. These are summarized in Table 4. The first routine allows for directly implementing σv θ -especially useful if the angular average can be performed analytically. By setting the Boolean variable

Routine (substitutes) Implements Returns
GetAnnMatrixAnalytic σv θ σv θ (PrepANNtheta) (as matrix) GetSecondMomentScat W Ω C el 2, (PrepSCATT for cBE) C el 2 Lastly, if the Boolean variable KDonly is true, annihilation is switched off, and one can investigate kinetic decoupling only (see, e.g., Ref. [37]) with the cBE and fBE routines. All preparing routine calls for annihilation are then redundant. For the fBE routine, this also allows to check accuracy settings, see Section 3.3.

Accuracy control parameters
The accuracy of preparatory and Boltzmann solver routines is controlled by a set of global parameters. In particular, the accuracy parameters introduced in Section 3.2 translate to DRAKE internal names as follows. The local maximum error, err i , is stored in nBEerr, cBEerr and fBEerr. Code-internal names for errNewton (n+1) i (controlling the Newton iterations) are cBEerrNewton and fBEerrNewton. The maximum relative change for any Newton iteration depth (n), errMaxNewton i , is set by cBEerrMaxNewton and fBEerrMaxNewton. The default values of these variables are generally not expected to require adjustments.
The number of phase-space density elements used by the fBE routine is set by the integer variable qN. For dimensionless momenta q = p/ m χ T in the default interval between qmin = 10 −8 and qmax = 12, qN ranges between about 120 and 400 in the pre-implemented models. Note that the runtime of fBE scales quadratically with qN, implying that it is worth optimizing the choice of this parameter in a newly added model. This can be done by, e.g., checking for numerical artifacts of too low phase-space density resolution, as described in Section 3.3.
We turn now to the accuracy parameters for the preparatory routines. PrepANN and PrepANN2 compute σv and σv 2 , respectively. To do so, these routines adaptively tabulate (ln x, ln σv (2) /GeV −2 ) pairs in a pre-set time (x) interval. The initial, minimum number of table entries in this time interval is governed by the parameter Nx. The adaptive procedure for increasing the tabulation density is based on the comparison of middle values, obtained from first and third order interpolations of each two consecutively tabulated points. The variable iacc and imax control this adaptive tabulation by setting the maximum allowed relative difference between those middle values and the maximal number of refinement iterations allowed, respectively. pg sets the precision goal (in digits) for the intrinsic Mathematica function NIntegrate that is used for numerical integration. After this construction, the tables are first interpolated and then exponentiated, resulting in the functions svx[x] and sv2x[x] listed in Table 3. Wolfram Language adaptive mesh procedure (used in plotting routines). The square root of the initial number of grid points is plotPoints. The recursion limit for grid adaption is controlled via plotMaxR. This irregular grid is then converted, by interpolation, into the regular, much denser grid tpregsvtheta that is used in the actual calculations. The square root of the number of points in this regular grid is set by Np, which can require values up to several thousand to maintain sharp features. The conversion time is rather short (compared to the overall fBE runtime), though scaling quadratically with Np. 13 Effective degrees of freedom The present release of DRAKE uses the SM effective number of energy and entropy degrees of freedom from Drees et al. [79] in form of a table located in the source directory. Alternative tables (in the same format) can be used by adding a corresponding file to the source directory, and replacing the string in dof="dof_Drees_etal.dat" in <settings> with this new file name. For other degree of freedom calculations, see the recent Ref. [80] and references therein. 13 fBE needs for every time step a qN 2 -sized matrix, whose entries are σv θ as a function of q i and q j . This matrix is provided by the sub-routine GetAnnMatrixTab, which linearly interpolates the regular tpregsvtheta in natural logarithm space. The regularity ensures fast access time, which overall saves much more computational time than needed for the conversion of the irregular to this regular table.

B: Elastic scattering beyond Fokker-Planck approximation
The DRAKE routines cBE and fBE can compute the relic abundance with the elastic scattering collision term C el once the model dependent quantity W Ω ≡ 1 4π dΩW , with W as introduced in Eq. (8), is provided. To illustrate this, we discuss concrete examples in Section B.1 where many of the integrals appearing in W Ω can be analytically performed, resulting in fairly simple expressions. For those cases, we also make a direct comparison between the resulting C el and its Fokker-Planck approximation. Finally, in Section B.2 we explain how C el is implemented in the DRAKE code and how to adjust it to any given W Ω .

B.1: Examples
As a first example, we consider non-relativistic heat bath particles, with g ± (ω) [ e −[m f +k 2 /(2m f )]/T , and scattering amplitudes independent of Mandelstam variable s. In this case most of the integrals entering in W Ω can be solved analytically, leading to We note that this expression is valid also for relativistic DM, and for an arbitrary dependence of |M| 2 χf ↔χf on t = q 2 0 − q 2 , where q = |p − p| is the momentum transfer and q 0 =Ẽ − E the energy transfer. The remaining integral can be performed analytically for amplitudes with a simple power law dependence on the Mandelstam variable t. In particular, in the case of a constant scattering amplitude, as in the sub-threshold annihilation scenario discussed in Section 4.3, this can be further simplified to: with erfc(x) being the complementary error function and As our second example, we consider again a scattering amplitude independent of s, but this time for ultrarelativistic heat bath particles, with g ± (ω) [1 ∓ g ± (ω)] e −k/T (1 ∓ e −k/T ). In that case, W Ω can be reduced to two remaining integrals: As before, this holds also for relativistic DM, and an arbitrary t-dependence of the scattering amplitude. In the above expression, the energies are given by ω = We now turn to explicit comparisons between C el [f χ ] calculated in these limits and the Fokker-Planck approximation C FP [f χ ]. For definiteness, and to have a non-vanishing scattering term, we choose a classical DM phase-space distribution f χ = e −E/Tχ , but with T χ slightly different from the heat bath temperature T . We will further always consider the product of p 4 /E 2 with the collision terms in these comparisons, because it is dp(p 4 /E 2 )C el,FP that determines the scattering strength C el 2 in Eq. (11), giving the y(x) evolution in the cBE approach.
We start by comparing, in Fig. 9, C el with W Ω from Eq. (B.2) to C FP with γ(x) from Eq. (31). The particular combinations of r (= m f /m χ ) and x values that we use here, for illustration, are chosen such that i) rx 1 to make sense of the non-relativistic bath assumption, and ii) γ/H ∼ 1 (close to kinetic decoupling) for |M| 2 χf ↔χf = λ 2 with λ = 0.05. At late times, i.e. for large x values and hence r 1, both scattering collision term descriptions are very close to each other (right panel). This is expected because small momentum transfer is guaranteed in this regime. For r 1 on the other hand, relevant for sub-threshold annihilations, small momentum transfer is not guaranteed. For example, for r = 1.1 (left panel) the Fokker-Planck approximation underestimates the collision term up to about 30 % at x=20 compared to the more accurate C el description (at these particular values of y eq /y). Similar differences are shown in Fig. 10, where instead the integrated quantities dp(p 4 /E 2 )C el,FP , as entering cBE, are compared. This qualitatively explains the differences in the relic-abundance results presented in Section 4.3: C el is more efficient compared to C FP in keeping DM in kinetic equilibrium, thus leading to a relic density closer to that of the nBE approach. C el non-relativistic bath C FP non-relativistic bath y eq /y=1.2 y eq /y=0.8  )C FP dp / ∫ (p 4 /E 2 )C el dp Fig. 10: The magnitude of the second moment of the Fokker-Plank approximation relative to the full scattering collision term, as a function of the DM temperature, for a constant scattering amplitude and non-relativistic bath particles. We note that for thermally produced DM y eq /y < 1 can only be achieved through (annihilation) processes that lead to DM self-heating.

20
to the corresponding Fokker-Planck approximation in Fig. 11 and Fig. 12, respectively. For low temperatures (right panels), as expected, the Fokker-Planck scattering term is always an excellent approximation. In the hightemperature regime (left panels), on the other hand, less separated scales (the momentum transfer scale ∼ T versus typical DM momenta ∼ m χ T χ ) can become an issue for the Fokker-Planck approximation -especially if amplitudes support larger momentum transfers as, e.g., in the |M| 2 χf ↔χf = −t/Λ 2 example. We conclude that the error induced when adopting the much simpler Fokker-Planck approximation, at temperatures relevant for chemical freeze-out, depends on the exact form of the scattering amplitude. The potential impact on the relic density in scenarios with early kinetic decoupling is therefore also unavoidably modeldependent. However, as for example in the threshold scenario discussed in Section 4.3, we expect that the dominant effect of early kinetic decoupling on the relic density can in many cases still be fairly well captured by the Fokker-Planck approximation.

B.2: Implementation
In discretized form, the full elastic collision term can be written as a linear operator of the form C el [f χ ]/E = C el · f with f = (f 0 , ..., f N ) T being the momentumdiscretized phase-space density f χ . The matrixĈ el can be written on a momentum grid, with p,p ∈ {p 0 , ..., p N } and equal spacing ∆p, as: GetScatMatrix and GetSecondMomentScat, respectively (see also Table 4). For the sub-threshold annihilation case, these two routines can be found in the model file TH.wl, where Eq. (B.5) with W Ω as in Eq. (B.2) is implemented. To switch from the Fokker-Planck approximation to full C el in the cBE and fBE approaches, set FullCel=True in the setting file. These two routines can be adopted to other scattering scenarios, by replacing the W ij dependent parts inside the routines.