Impact of dark matter self-scattering on its relic abundance

Elastic self-scatterings do not change the number of dark matter particles and as such have been neglected in the calculation of its relic abundance. In this work we highlight the scenarios where the presence of self-scatterings has a significant impact on the effectiveness of annihilation processes through the modification of dark matter momentum distribution. We study a few example freeze-out scenarios involving resonant and sub-threshold annihilations, as well as a model with an additional source of dark matter particles from the decays of a heavier mediator state. Interestingly, when the calculation is performed at the level of dark matter momentum distribution function, we find that the injection of additional energetic dark matter particles onto the thermal population can lead to a $\textit{decrease}$ of its final relic abundance.

Elastic self-scatterings do not change the number of dark matter particles and as such have been neglected in the calculation of its relic abundance. In this work we highlight the scenarios where the presence of self-scatterings has a significant impact on the effectiveness of annihilation processes through the modification of dark matter momentum distribution. We study a few example freeze-out scenarios involving resonant and sub-threshold annihilations, as well as a model with an additional source of dark matter particles from the decays of a heavier mediator state. Interestingly, when the calculation is performed at the level of dark matter momentum distribution function, we find that the injection of additional energetic dark matter particles onto the thermal population can lead to a decrease of its final relic abundance.

I. INTRODUCTION
The extremely dense and hot plasma of the very Early Universe exhibits exceptional conditions for particle creation. Every new particle species that i) is coupled to the Standard Model (SM) states, ii) has a non-vanishing mass not exceeding the reheating temperature and iii) is stable on cosmological timescales, will inevitably end up with a relic thermal population contributing to the present-day dark matter (DM) density. In the simplest scenarios this thermal component accounts for all of the observed DM, with Ωh 2 = 0.120 ± 0.001 [1]. If adopted, this requirement puts significant constraints on the parameters of a given DM model that affect the rates of the particlenumber-changing processes.
An attractive, yet not overly restrictive assumption is that the interactions between the DM particles and SM plasma are sufficiently frequent to enforce chemical equilibrium at some time in the very Early Universe. In such scenario the relic density of DM is predominantly determined by just one quantitythe cross section of annihilation -in a wide range of DM masses. Effectively, when the rate of annihilation drops below the rate of the Universe expansion, the dark sector departs from the chemical equilibrium with the SM plasma and the number of DM particles in the comoving volume ceases to change with time (freezes-out), hence establishing the relic population. The general approach to determine the relic abundance of DM is solving the Boltzmann equation (BE) that describes the evolution of the DM distribution function in the expanding Universe. To calculate the relic abundance often only the 0th moment of this equation, tracing only the particle-number density, is considered, which is what is used in numerous existing numerical packages, e.g. [2][3][4].
In a broader class of DM models thermodynamics in the early Universe can be more elaborate and the relic abundance can be generated in other ways than described above, but the leading principle that the relic density is determined by the interplay of DM number-changing processes remains. * andrzej.hryczuk@ncbj.gov.pl † maxim.laletin@ncbj.gov.pl However, the rate of number-changing processes depends not only on the interaction strength, but also on the characteristics of the DM population, in particular its number density and momentum distribution. While the former is determined by chemical equilibration and decoupling governed by annihilations, the latter is related to the local kinetic equilibrium that is maintained mostly by elastic and inelastic scatterings on the particles from the SM plasma. Although in the typical models of weakly interacting massive particles (WIMPs) the kinetic equilibrium is maintained long after the freeze out, exceptions to this standard scenario exist in even simple models [5][6][7][8][9][10][11][12] and are expected to occur much more often in more involved scenarios containing processes actively disrupting local thermal equilibrium, e.g. decays of heavier states or self-heating [13] due to semi-annihilations, cannibalization [14] or conversions (see e.g. [8,15,16]). In such cases the crucial question is whether the scattering processes are efficient enough to force the momentum distribution to follow the equilibrium one with the same temperature as the SM plasma or not.
The impact of elastic scatterings between DM and the thermal-bath particles has been studied in recent years both at the level of tracing the evolution of DM temperature, see e.g. [5][6][7][17][18][19][20], alongside the number density in a coupled systems of Boltzmann equations (cBE) for the 0th and the 2nd moments of the distribution function, as well as at the level of numerical solution of the full momentum-dependent Boltzmann equation (fBE) [7,12,21,22]. When comparing these two approaches it was noted [7] that although cBE and fBE give in general consistent results, they can substantially differ in scenarios with strong velocity-dependent annihilations, e.g. due to a resonance or a threshold. This is because cBE enforces the shape of the distribution to be of the Maxwell-Boltzmann form, which is equivalent to insisting that DM self-scatterings, redistributing energy in the dark matter component, are very efficient. While without the explicit inclusion of DM self-scattering on top of annihilation and elastic scatterings on the thermalbath particles, the fBE method effectively neglects all such processes altogether. Hence, which of these approaches gives a better estimate clearly depends on the actual strength of self-scattering processes.
In this paper we implement, for the first time in the literature, the complete DM self-scatterings at the level of fBE in the context of thermal relic density calculation 1 and investigate how do they modify the energy distribution of DM and ultimately its relic density for three different models of DM with strong velocity dependence of annihilation processes. We compare the results of the different approaches discussed above to shed light on their ranges of applicability regarding the rate of self-scattering. In particular, we explore a model with the injection of an energetic component of DM in the form of heavier particle decay products.
It is worth mentioning, that as far as elastic scatterings on thermal-bath particles and annihilations are typically strongly tied together in the underlying DM model, the self-scattering processes are often unrelated and in particular can have much larger cross sections. In fact, models of strongly interacting DM have gained a substantial attention in the literature as a possible solution of the discrepancies between observations and theoretical prediction for the density profiles of DM in small scale structures (see e.g. Ref. [24] for a review).
The paper is organized as follows. In Sec. II we describe our implementation of self-scattering processes in the Boltzmann equation. Sec. III discusses applications to two models with the standard freeze-out production, while in Sec. IV we extend the analysis to an example with an additional non-thermally produced component and its interplay with the freeze-out mechanism. Finally, Sec. V concludes.

II. DARK MATTER SELF-SCATTERING
Before we discuss the specifics of self-scattering processes, let us briefly review the formalism used in the calculations of the thermal relic density. The evolution of the DM component is well described in the semi-classical limit by the Boltzmann equation that in Friedmann-Roberston-Walker space-time takes the following form: where f χ (t, p) is the DM distribution function depending on time t (in what follows replaced with temperature T ) and momentum p (or equivalently, energy E), H stands for the Hubble expansion rate and C denotes different collision terms that are relevant in the Early Universe. As an example we provide here the structure of the collision term for a particle i that 1 For an implementation within a relaxation-like approximation see Ref. [23].
participates in a general 2-to-2 process where g i is the number of degrees of freedom of particle i, ω andω (4-momenta k andk) are the energies of the final state particles m and n, and E andẼ are the energies of the initial state particles i and j respectively (4-momenta p andp). The amplitude squared |M| 2 is summed over both initial and final internal degrees of freedom. The corresponding distribution functions are marked with the respective indices. The first term in this expression is generally referred to as the gain term, while the second -the loss term. The signs in front of the distribution functions of the products in the forward and backward reactions depend on the spin statistics describing these particles. If these particles constitute a very dilute gas [1 ± f ] ≈ 1 and these factors can be neglected.
The standard approach [25] of solving only the single Boltzmann equation for the number density n χ (nBE) can be obtained from Eq. (1) by the integration over the momentum p leading to (3) The collision terms for elastic and self-scattering do not change the number density and therefore cancel out after the integration. However, in order to solve this equation one needs to know the form of f χ (t, p), so that the r.h.s. of the equation can be integrated and expressed in terms of the number density. A common assumption is that the distribution of DM has an equilibrium shape (Fermi-Dirac/Bose-Einstein distributions in general or Maxwell-Boltzmann distribution in the dilute or non-relativistic limit) that corresponds to the temperature of the SM plasma and with a potentially non-zero chemical potential that is effectively solved for. This assumption is often justified since the elastic scattering processes on SM particles typically proceed at a large enough rate. In cases when the elastic scatterings cannot maintain local thermal equilibrium, but the shape of the distribution function is still close to the thermal one, albeit with T χ = T , the cBE approach is expected to give an accurate prediction not only for the DM relic abundance, but also for its temperature evolution. This system of equations for the number density and temperature is obtained from Eq. (1) by the integration over (g χ /(2π) 3 ) d 3 p/E and (g χ /(2π) 3 ) d 3 p p 2 /E 2 respectively. In the following we show the results obtained with these two approaches only as a comparison to the full treatment we study in this work. Thus, for more technical details regarding nBE and cBE we refer to Ref. [7], while below we discuss our implementation of the fBE, and especially the self-scatterings.
The collision term for annihilation of DM particles into two SM states (neglecting the [1 ± f χ ] factors) is given by where g stands for the thermal distribution for a considered thermal-bath state. The gain term in the annihilation collision term does not contain any unknown distribution functions and in principle can be calculated explicitly. The loss term contains two DM distribution functions: one of them can be taken out of the integration, the other has to be integrated over the corresponding momentum, while the residual part of the integrand can be expressed in terms of annihilation cross section. To perform this integration numerically the unknown distribution function can be regarded as a combination of discrete components f χ (p i ) for a given value of momentum p i . Thus, the momentum-dependent Boltzmann equation is split into a system of ordinary differential equa-tions for each momentum component and the integration is approximated with a weighted sum of these components. This scheme in particular is realized in the DRAKE code [12], which we use for the solution of the BE, except for the self-scatterings (see below) which are not included in the current public version of DRAKE. The same procedure can be in principle applied to the loss term of the elastic collision term which has the same structure as the gain term, but with the energies transformed as indicated. However, in the limit of small momentum transfer the whole C el can be expressed through the DM distribution function and its derivatives without any numerical integrations left [26] (see Appendix B). Finally, the decay collision term C dec can be simplified to an analytical expression (see Sec. IV), as long as the decaying particle is described by an equilibrium distribution. The collision term for self-scattering has the following general expression (neglecting the [1 ± f χ ] factors) The first term (both loss and gain parts) describes the scattering on particles and the second term the scattering on antiparticles. The factor 1/2 in front of the first one takes into account the symmetry between the identical particles χ with momenta that are integrated over. In the absence of CP -violating processes in the dark sector the distribution function for particles and antiparticles is always the same, thus the self-scattering collision term can be written in terms of an effective amplitude squared and one set of gain and loss terms. In comparison to elastic scattering, the rate of selfscattering is suppressed by an additional f χ in the collision term, especially with respect to light SM states that are greatly more abundant in equilibrium. However, it is an insufficient reason to claim that selfscatterings play a little role in shaping the DM energy distribution. First of all, an average relative momentum transfer for elastic scattering is δp/p ∼ (T /m χ ) 1/2 1, while for self-scattering δp/p ∼ 1, so an effective energy redistribution in the latter case does not require many collisions. It is useful to com-pare the characteristic relaxation times for both processes τ r ∼ N coll /Γ, where Γ = n χ σv is the rate of scatterings and N coll is the number of collisions required to substantially change the momentum of DM (see also Ref. [26]). In case of elastic scatterings N el coll ∼ m χ /T , the density of relativistic SM particles n SM ∝ T 3 and the overall relaxation time scales with temperature as T −6 [27]. For self-scatterings after freeze-out the density of DM particles scales with temperature by the same law due to the expansion, but N self coll ∼ 1 and for the models that we consider in Sec. III A and IV σv self ∝ T in the non-relativistic limit. Thus, the relaxation time for self-scattering scales with temperature as T −4 , which means that self-scattering processes remain an effective mean of equilibration longer than elastic scatterings.
Secondly, self-scattering can rely on different couplings (or even on a different type of interaction) than the ones that govern elastic scattering. For instance, in the case of scalar DM model considered in Sec. III B self-interaction can arise from a simple φ 4 vertex interaction with the amplitude squared proportional to the square of the respective coupling, while elastic scattering on SM fermions is loop-suppressed. In the case of fermion DM coupled to a vector mediator (Sec. III A), the rate of elastic scattering can be suppressed by the squared ratio of the two couplings, given that the coupling of the mediator to the heat bath fermions is smaller. In other models of DM self-scattering can be boosted w.r.t. elastic scattering by an s-channel resonance, Sommerfeld enhancement, etc. In addition, DM self-interactions are generally not as constrained by observations as the elastic scatterings. For example, the upper bound on the cross section of electron scattering for a DM particle with the mass of 1 GeV is ∼ 10 −34 cm 2 [28], while the cross section of DM selfscattering for that mass can be as large as ∼ 10 −24 cm 2 [29].
From a technical point of view, the additional complication introduced by the self-scattering collision term in Eq. 6 comes from the gain term that contains two distribution functions of the products. While the loss term can be treated in the same way as for annihilation, the gain term for self-scattering cannot be simply formulated in terms of the cross section and the presence of two unknown function in the integrand leads to complicated angular dependencies during the integration (see Appendix A). Since selfscatterings are not implemented in the current version of DRAKE, we merged the existing version with the program for the calculation of the self-scattering collision term that we developed. Though the loss term in Eq. 6 can be expressed through the self-scattering cross section and implemented numerically with one summation of the momentum components, we use the same procedure as for the gain term to achieve a better numerical cancellation between the two terms close to equilibrium point, since the same interpolation procedure is used in both cases.

III. IMPACT IN FREEZE-OUT MODELS
During the freeze-out process the self-scattering does not introduce any direct change in the DM number density, but indirectly it can significantly affect the annihilation rates. To exemplify this we have chosen to present results of a study of two models introduced in Ref. [12]: the generic vector resonance and sub-threshold scenarios. In both cases the kinetic decoupling and non-equilibrium shape of f χ (p) can have a strong impact on the final relic abundance and as noted in Ref. [12] there can be a substantial difference between the cBE and fBE approaches. As it was pointed out in Ref. [7] it is expected that inclusion of self-scatterings to fBE treatment should in the limit of large self-interactions lead to result coinciding with cBE. In this section we demonstrate that this is indeed the case and quantify how strong the self interactions need to be to have an impact.
A. Vector resonance model The arguably most common scenario where the DM annihilation cross-section has a strong velocity dependence arises in models with s-channel resonance. For concreteness let us take exactly the same model as in Ref. [12] where the resonance is mediated through an exchange of a generic vector mediator A µ , with the interaction Lagrangian The model can be described by a set of five parameters: the DM mass m χ , the mediator mass m A , the mass ratio of heat-bath fermions to DM r ≡ m f /m χ and finally the coupling constants λ f and λ χ . Out of these input parameters it is convenient to define deviation from the exact resonance position δ ≡ (2m χ /m A ) 2 − 1 and a dimensionless measure of the total decay width of A µ ,γ ≡ Γ A /m A . Note, that compared to the discussion in Ref. [12] we separate the couplings λ f and λ χ , as the self-interactions break the degeneracy between them in the calculation of the DM annihilation cross section. Indeed, the annihilation cross-section for the process χχ → A → ff , can be written as [12] s is the center-of-mass energy, and being the Breit-Wigner propagator. 2 The selfscattering amplitude squared in this model consists of two contributions, as in Eq. (7), but both of them depend only on the coupling λ χ wheret ≡ t/(4m 2 χ ),ũ = 1 −s −t and β 1 and β 2 are functions ofs andt defined in the Appendix B.
In the calculations we also include elastic scatterings on the thermal-plasma fermions f , which is exactly in the same way as in Ref. [12], to where we refer the reader for more details.
Let us start the discussion of the results from presenting the evolution of the yield Y = n χ /s and temperature parameter y = m χ T χ s −2/3 , where s(T ) is the entropy density, for a benchmark point with the resonance having SM Higgs-like widthγ = 3 × 10 −5 , and relatively heavy annihilation products, r = 0.5 with m χ = 100 GeV. Fixing the widthγ and requiring that the nBE solution provides the observed relic abundance defines the couplings λ f = 10 −3 and λ χ = 5.85 × 10 −2 and ultimately the strength of selfscattering process. The result of the evolution is given in Fig. 1 for the nBE (blue), cBE (green), fBE without (orange) and including self-scatterings (black). For this benchmark δ = −0.05 and thus the resonant annihilations deplete momenta around the peak of the distribution at the time of freeze-out, resulting in a slight temperature raise at first and then an abrupt chemical and kinetic decoupling. Self-scatterings reshuffle the DM particles' momenta, re-populate the regions depleted by annihilation and thus prolong the freezeout, leading to a lower final abundance. This can be directly seen in Fig. 2, where the four time snapshots of the momentum distribution are given. Comparing to the result obtained with self-scatterings (black), which make f χ (p) retain shape close to thermal, the curves without self-scattering (orange) show a significant dip in the distribution for the momenta that are slightly above the peak of the distribution. The particles with these momenta are efficiently depleted by the resonant annihilation, while elastic scatterings on the thermal bath are not sufficient to replenish them effectively. The overall distribution is also visibly shifted from the equilibrium one (blue) indicating that T χ is significantly lower than T indeed. The bottom panel of Fig. 2 highlights the relative size of the difference between the two distributions, with different colours signifying different time snapshots. One can see that the self-thermalization due to selfscatterings in this model is rather efficient, in large part due to the fact that the distortion introduced by v-dependent annihilation is limited. Nevertheless, even such relatively small deviation from thermal shape visibly affects the relic abundance through modifying the annihilation rate. Quantitatively, Fig. 3 shows the relative difference between the thermally averaged cross sections σv calculated with the resulting f χ (p) and with a thermal distribution of the same temperature. The result without self-scatterings (orange line) around the freeze-out time x ∼ 25 shows significant deviation from the thermal one, while with self-scatterings (black line) a much milder one. When the temperature drops, to about x ∼ 50 and further, both solutions predict an enhanced annihilation rate, but with somewhat different dependence. In particular the peak of this enhancement happens earlier for the solution with self-scatterings than without, which is a consequence of an impact on resonant annihilation by an interplay between the temperature drop and the distribution shape modification. The deviations from unity thus show that even with the self-scatterings the annihilations introduce too much of a disruption to allow maintaining an equilibrium shape. Note that these are both normalized to T χ corresponding to the given fBE solution in order to highlight the effect of the non-thermal shape only. This can be contrasted with the blue line including the impact of the temperature change as well, where the comparison is made to σv T .
Finally, Fig. 4 shows how the change of the coupling λ χ affects the relic abundance for the same benchmark scenario. If this coupling is not fixed by the relic abundance requirement, then fixingγ introduces relation λ f (λ χ ) and one can vary the strength of self-scatterings. Note though, that there is a maximal value, in case of this benchmark point λ max χ = 5.855 × 10 −2 , above which it is not possible to obtain widthγ = 3 × 10 −5 . In order to highlight the most interesting region, in which self-scatterings are as effective as possible, the x-axis of the figure displays the distance from this maximal value λ max χ . Top panel shows the result for the relic density, while the bottom one -the ratios of the result for Ωh 2 obtained with fBE without (orange) and with self-scatterings (black) to the cBE one. For small values of λ χ (on the right edge of the plot) the fBE result coincides with the one without self-scatterings whatsoever, while for the val- ues approaching the maximal one (on the left edge) it departs towards the result of cBE. All in all, these results suggest that for typical values of the self-scattering strength the fBE as currently implemented in DRAKE gives a better approximation of the actual result than the cBE approach. However, this statement is model dependent and when a precise result for the relic abundance is called for, one should in principle fully include self-interactions.

B. Sub-threshold model
A second common scenario where the DM annihilation cross-section has a strong velocity dependence is when the DM annihilation process has a threshold at s > 4m 2 DM . Below this threshold the annihilation can be kinematically impossible, sometimes dubbed 'forbidden' DM [30], or non-zero, but suppressed. Again for concreteness let us take exactly the same model as the sub-threshold model in Ref. [12], composed of two scalar particles, where φ 1 takes the role of the DM, while φ 2 is in thermal contact with the heatbath fermions f . The interaction Lagrangian is given by where compared to the discussion of this model in Ref. [12] we added the self-interaction term, that was not implemented previously. We will assume that the scalars are close in mass and the DM is slightly lighter, i.e., r ≡ m 2 /m 1 1. In such regime, to the lowest order, the total DM annihilation takes place through the process φ 1 φ 1 → φ 2 φ 2 and the cross section is given by σv lab = g 2 32π while the momentum transfer rate between the DM and the heat-bath fermions f is strongly suppressed [12]. The amplitude squared of self-scattering in this case is simply a constant and equal to λ 2 .
In the resonance example we discussed in detail one benchmark model point -here instead let us focus on the relic density as a function of self-scattering strength for a representative set of parameter points. In Fig. 5 we show the effect on the relic density for four values of r ∈ {1.001, 1.1, 1.15, 1.2} using cBE (dotted), fBE without self interactions (dashed) and full calculation (solid) as a function of the self-coupling λ. Increasing the self-scattering strength makes the full result go from being the same as fBE λ = 0 to approaching the cBE result, which agrees with the expectations. However, performing the analysis at the quantitative level reveals that, at least for a forbiddenlike sub-threshold model at hand, one would require λ > 1 to actually significantly depart from the fBE λ = 0 result. This confirms the observation from the previous section, that the fBE result from DRAKE is expected to be typically a better estimate for the relic density, than the cBE one.
Before closing this section a comment on the technical side of the numerical computation is in order. Eq. (6) after discretization takes the form of a 3 di- mensional matrix of entries being 2 dimensional angular integrals, see Eq. (A1). Therefore, the mitigation of numerical inaccuracies by discretizing on a denser grid is rather costly CPU-wise with O(N 3 ) scaling with the grid density for tabulation and also O(N 2 ) scaling for generating the collision term matrix at every x-step of fBE. Even though typically this matrix is rather well-behaved, this problem becomes especially relevant when the annihilation predominantly relies on the high momentum tail of the distribution. This is the case for the sub-threshold model where the tail is much more prone to numerical error than the bulk of the distribution. All in all, the results presented on the plot in Fig. 5 required several hours of CPU time per point and still have small, but visible irregularities in the full result at large values of λ. This attests to the level of numerical accuracy that one can achieve within a manageable CPU cost with our current implementation. 3 3 The numerical code we used is available on request and is planned to be publicly released as an additional optional package to a future DRAKE version.

IV. SELF-THERMALIZATION OF A NON-THERMAL COMPONENT
The effects of self-scattering are expected to be particularly important for scenarios in which a considerable portion of DM is produced non-thermally on top of the thermal component. Quite commonly particle physics models with the DM candidate(s) contain heavier particles that can decay with a production of one or several DM states. If the lifetime of these heavy particles is sufficiently long, the contribution from the decays does not simply annihilate away and return the distribution to the equilibrium, but can noticeably alter the evolution of DM distribution, its density and other properties [17,[31][32][33][34]. This additional non-thermal contribution can also arise, for example, from the bubble collisions following a first-order phase transition [35] or primordial black holes evaporation [36]. If the velocity-averaged cross section of DM annihilation is essentially momentum-independent this injection will have an impact on the rate of annihilation solely by the increase of the density -the resulting relic density will be determined by the interplay of the prolongation of annihilation that depletes the density and the continuous supply of new particles that increases the relic abundance. However, if the annihilation is strongly velocity dependent the effect of the distribution on the annihilation rate is more complicated. If the injected component is rather energetic w.r.t. the thermal one, self-scattering processes will lead to the redistribution of DM particles into the region of the phase space with a larger momentum and hence it can noticeably affect the velocity-averaged cross section. Moreover, the effects of self-scattering on the energy distribution of DM with a non-thermal component can have consequences that go beyond just the impact on the relic density. For instance, the shape of the relic distribution of light DM particles at later stages of the evolution of the Universe can affect the formation of large scale structure (e.g. [37][38][39]). However it is worth exploring, the analysis of this phenomenon stays out of the scope of our paper.

A. Model setup
To study the self-thermalization of a non-thermal component we take a sterile-neutrino-like model of DM that is coupled to a scalar singlet field S, which has been previously studied in the literature (e.g. [38,40,41]), with an additional U (1) gauge interaction (dark electromagnetism). In Ref. [22] a similar model is considered in the context of the impact of non-thermal processes on the DM distribution function. However, it does not include the additional U (1) making the self-scattering processes absent in their case. The Lagrangian of the model looks as follows where is the mixing parameter between the photon and the dark photon A µ , D µ = ∂ µ −ie A µ and V (S, H) is the Z 2 symmetric part of the scalar potential The reason we impose the nearly exact, as explicitly broken only by the interaction of S with DM fermions, S → −S symmetry is twofold. For one it motivates the choice of the effective Yukawa coupling y to be small and lead to a long lifetime of S (we will focus on a case where m S is sufficiently large to allow the decays to a pair of energetic DM particles). But also it allows for a more clear discussion of the process studied in this work without being sidetracked to other, well known effects. In particular, significant explicit or spontaneous breaking of this Z 2 symmetry would lead to terms in the Lagrangian that would allow for (loop suppressed, though typically quite efficient) S decay to the SM states. This would not affect directly the discussion of self-thermalization and the size of the DM component coming from S decay could always be adjusted by modifying the S freeze-out. Nevertheless, it could lead to a substantial entropy increase due to the decays of S to radiation and introduce another impact on the relic density of DM that has been already very well studied in the literature (see e.g. [42][43][44]).
In the case that we consider the mass of the dark fermion m χ is introduced explicitly by the mass term. We do not stipulate here the exact mechanism by which the dark photon gets its mass or a UV complete theory of the mixing between the two gauge interactions, however several possible scenarios have been studied in the literature and we refer the reader to a recent review of these models [45].
We focus on a region of the model parameters that satisfy the following conditions: a) the coupling of the Higgs field to the singlet scalar S is strong enough to keep it in kinetic equilibrium with the SM bath around the DM freeze-out; b) m A < 2m χ , so that the dark photon can only decay to SM states. Then the contribution to the collision term of χ(p) from the decay of S(k) is given by (16) where we have neglected the inverse decay process, as is appropriate for a long-lived S. Using the fact that f S (ω) ∝ f eq S (ω) and that f χ (E) 1 the term above simplifies significantly and is completely independent on the DM distribution function. For the decay amplitude squared after performing the integrations one arrives at The normalization of this decay term is proportional to the number density of S particles, which can be in or out of equilibrium. It is in turn determined by the chemical decoupling of S from the thermal bath that is governed by its annihilation processes. In our numerical implementation we solve an nBE-type Boltzmann equation for its evolution including the decay process, but neglecting back reaction of χ, which is a very good approximation as long as the decay happens somewhat later than its decoupling. The self-scatterings and elastic scatterings are mediated by the dark photon A µ in the same fashion as in the model considered in Sec. III A. Depending on the ratio of the masses of DM and the dark photon m χ /m A this model can reproduce both of the DM annihilation patterns considered in Section III. In the region of masses m A /m χ 2 the DM annihilates to SM electrically charged states via the resonantly-enhanced s-channel mixing between the two photon mediators (resonance regime). In the region m A /m χ 1 annihilation channel to two dark photons is opened when s > 4m 2 A . Far from any resonance, annihilation cross section to dark photons is proportional to e 4 , while annihilation to SM states via photon mixing is suppressed by the factor 2 α/e 2 , where α is the fine structure constant. Since the upper limit on the kinetic mixing for the value of the dark photon mass that is relevant for our study is ∼ 10 −3 (see e.g. Fig. 3.3 in Ref. [45]), the difference of annihilation rates below and above the threshold of s = 4m 2 A is significant, effectively leading to the sub-threshold regime. Despite the smallness of the dark photon remains in chemical equilibrium throughout the freeze-out of DM.
Finally, before discussing the results a comment on the elastic scatterings on bath particles is in order. We follow the implementation in DRAKE, briefly summarized in the Appendix A, which was derived in a semi-relativistic Fokker-Planck approximation applicable to the thermal freeze-out. In the presence of an additional relativistic component formally this treatment breaks down and the implemented elasticscattering term may deviate from the actual one. Nevertheless, as we will see, for the results shown the momentum transfer rate in elastic scatterings, γ(T ) is subleading compared to the other processes and therefore the current implementation is expected to be sufficient.

B. Results
Below we show numerical results only for the subthreshold regime, because the resonance regime for the given model with the constraints on requires an extremely sharp resonance to not overproduce the relic density, and so it is phenomenologically less interesting within the model at hand, while retaining similar qualitative behaviour. An example of the particle yield Y and the temperature parameter y is shown in Fig. 6 for a benchmark set of parameters m χ = 100 GeV, m A = 108 GeV, e = 1 and = 0.001, which is chosen such that the nBE approach with the decays switched-off (gray curve) reproduces the observed relic density. All the curves, except for the gray one, display the same inflection point after the density decouples from the equilibrium value. At this point the rate of particle loss due to annihilation and the rate of particle gain due to the decay are comparable. From this point the nBE curve (blue) grows somewhat, but the fBE curves (orange and black) decrease even further. While the full fBE (black), fBE without self-scatterings (orange), nBE at Tχ and SM plasma T (green and blue, respectively). Dotted blue and black lines show the momentum transfer and self-scattering rates respectively. Note that these are the rates are obtained using the actual solutions for nχ, not its equilibrium value. The raise of the annihilation rates in the fBE approaches starting at x ∼ 30 explains why the injection of extra DM particles can lead to a significantly enhanced annihilation and an ultimately decreased relic abundance. The green curve is shown to highlight the fact that it is not the change in temperature, but the shape of the distribution that is the main reason of this enhanced annihilation.
behaviour of the first curve is expected as the production of additional DM particles should increase the density, the fact that the account for the actual shape of DM momentum distribution in this case leads to the decrease of the DM density can seem quite surprising. However, it can be easily understood from the velocity-dependent annihilation pattern of the model and the momentum distributions that correspond to the different approaches. The evolution of f χ (p) is shown in Fig. 7. In the nBE approach the shape of the distribution is assumed to be unchanged even if the DM particles from S decay are in fact more energetic, hence the rate of annihilation is only slightly affected by the presence of additional DM particles and the decay gain term dominates over annihilation in the density evolution until it becomes too small to noticeably increase the abundance of DM. In the fBE case with the switched-off self-scatterings (orange) the decays create a small bump in the distribution function with a characteristic momentum that is sufficient to overcome the annihilation threshold and thus the rate of annihilations in the DM gas is significantly boosted. The effective temperature of DM rises w.r.t. the SM plasma temperature (orange line in the right plot in Fig. 6), however as the most energetic particles in the bump annihilate away the effective temperature drops below the equilibrium temperature due to the cooling caused by the expansion. At this stage the majority of annihilation processes happens below the threshold again, hence the thermally averaged cross section becomes very small and the DM freezes out.
In the presence of self-scattering the fBE curve (black) displays a steeper drop of the particle yield as the injected DM component heats up the DM gas via elastic collisions such that more DM particles have the energies to overcome the annihilation threshold. In the example that we consider the rate of self-scatterings is large enough compared to the rate of annihilations 4 in the absence of self-scattering (cf. Fig. 8) that the injected component can effectively transfer the energy to the thermal component before that additional component is annihilated and the energy stored in it is dumped into the SM plasma. This heat significantly increases the temperature of DM and keeps the annihilation rate larger than the Hubble rate for a longer time, so that the relic density is established later and gets a smaller value compared to the case when self-scatterings are switched off. The second inflection point in the abundance curve appears when the supply of DM particles from the decays essentially ceases and little heat is injected in the DM gas. In spite of that the rate of annihilation does not decrease as fast as in the absence of self-scattering, because these scatterings promote the clustering of particles in the phase-space region of higher momenta than the characteristic thermal scale and hence prolong the annihilation above the threshold. Note that the enhancement of annihilation that we consider here is only relevant for the early history of DM evolution -long before the effects of DM annihilation can leave an imprint on cosmological observables the DM gas is cooled down by the expansion of the Universe to an extent that annihilation proceeds only below the threshold.
To summarize, we have considered a practical example of a DM model with a late non-thermal component and demonstrated that the offset between the relic densities predicted by the nBE and fBE treatments can reach up to a few orders of magnitude and that the impact of self-scatterings on the result is crucial as well.
The above discussion also highlights that there is no good way of formulating a general and simple model independent criterion as to when the DM distribution's departure from equilibrium shape significantly affects the relic density. The final effect comes from an interplay of not only the expansion, annihilation, elastic and self-scattering processes, but also potentially other processes that disrupt the equilibrium, e.g. decays or annihilations of heavier particles into DM. Moreover, comparison of just the rates for these processes is insufficient, as they can be very efficient for some range of the momenta while not for others, as exemplified in the discussed models. A useful rule of thumb of when to expect a possibility of departure from kinetic equilibrium is if the rate of a process that disrupts kinetic equilibrium is larger than the rate of elastic scatterings between the DM and bath particles. As to the importance of careful inclusion of self-scattering processes one can expect it to be necessary when the result obtained with cBE and fBE methods differs to larger extent than the desired accuracy.
Before ending this section let us mention that although we have not studied phenomenological implications of the observed effects, other than the relic abundance, we did check that the presented benchmarks within the example models are feasible. In particular, the present day annihilation cross section is well below the current observational limits and A decay lifetime is short enough not to spoil the cosmic microwave background anisotropies, nor the Big Bang Nucleosynthesis.

V. CONCLUSIONS
In this paper we investigated the impact of the DM elastic self-scattering process on the evolution of its momentum distribution function and formation of the relic abundance. Building upon the DRAKE framework and code we implemented numerically the selfscattering collision term and applied it to two example models with thermal freeze-out and one model with an additional source of DM particles from decay of a heavier long-lived state.
We found that in all these cases the effect of selfscattering on the thermalization of the distribution function is large enough to bring visible changes in the effective annihilation rates and therefore final relic abundance of DM. In the freeze-out models our result does follow the expectation of interpolating between the relic abundance obtained using the coupled system of Boltzmann equations for number density and temperature (cBE) and the numerical solution of for f χ (p) without including the self-scattering collision term (fBE). This not only validates both of these approaches in their respective limits, but also shows that very large self-scattering rates are needed to recover the cBE result, thus suggesting, though in a modeldependent way, that the fBE is typically a more accurate approach.
In the case of an additional source of DM particles from a decay of a heavier state the disruption of the close-to-equilibrium shape of the thermal component that follows can be very significant, making self-scatterings a crucial ingredient in obtaining the accurate predictions for the relic abundance, temperature and the shape of the DM distribution function. Additionally, we found an intriguing feature that can only be uncovered when studying the evolution at the level of the distribution function: the injection of high energetic DM particles on top of the freeze-out thermal component can lead to a decrease of the resulting relic abundance. This effect arises in the scenarios with strong velocity-dependent annihilations and is significantly enhanced due to DM self-interactions. It would be interesting to consider phenomenological implications of such an effect in concrete DM models, which we leave for future work.
Though in this paper we focused on the relic abundance, the provided study also has consequences for the prediction of the matter power-spectrum for light DM candidates that have some non-thermal component. In such situations the precise shape of the distribution function is needed to accurately predict the size of density perturbations and as a consequence can affect, e.g., the warm DM mass limit [39]. Our work shows that self-scatterings can play a significant role in such calculations. and k 2 cm = s − (m χ − m f ) 2 s − (m χ + m f ) 2 /(4s) evaluated at s = m 2 χ + 2ωm χ + m 2 f .