{'fraction_non_alphanumeric': 0.037587961799447095, 'fraction_numerical': 0.01234606685096758, 'mean_word_length': 4.389857771757534, 'pattern_counts': {'":': 1, '<': 3, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Boltzmann hierarchies for self-interacting warm dark matter scenarios

We provide a general framework for self-interacting warm dark matter (WDM) in cosmological perturbations, by deriving from first principles a Boltzmann hierarchy which retains certain independence from a particular interaction Lagrangian. We consider elastic interactions among the massive particles, and obtain a hierarchy which is more general than the ones usually obtained for non-relativistic (as for cold DM) or for ultra-relativistic (as for neutrinos) approximations. The more general momentum-dependent kernel integrals in the Boltzmann collision terms, are explicitly calculated for different field-mediator models, including examples of a scalar field (either massive or massless) or a massive vector field. As an application, we study the evolution of the interaction rate per particle under the relaxation time approximation, and assess when a given self-interaction is relevant in comparison with the Hubble expansion rate. Our framework aims to be a useful starting point to evaluate DM self-interaction effects in the linear power spectrum, necessary to then study its evolution all the way to non-linear stages of structure formation, where certain DM interactions were proven to be relevant.


Introduction
In the standard cosmological paradigm, ΛCDM, Dark Matter (DM) has long been a necessary ingredient in the Big Bang model of the universe and in understanding its evolution since the early stages. While the evidence for its existence is implied by its gravitational effects in astrophysical, galactic and cosmological structures; understanding the nature and composition of this species is still an elusive subject [1][2][3][4]. Several attempts have been made to explain this phenomenon by macroscopic objects, yet a microscopic origin of the DM phenomenon by a new particle species remains as the most plausible hypothesis [3][4][5][6][7]. On the early stages of DM research, active neutrinos appeared as promising candidates for this particle species [8,9]. However, neutrinos are "hot" Dark Matter (HDM) with a free-streaming length which erases structures at large scales [10], while numerical simulations have shown that such "top-down" structure formation is incompatible with clustering constraints [11]. Cosmological data favored the adoption of the ΛCDM paradigm [12]: in the standard scenario, DM is assumed to be produced in a thermal distribution and modeled as collisionless after it decouples from the other species. The effective decoupling is assumed to occur at a temperature smaller than the DM mass so that the distribution corresponds to non-relativistic particles. The traditional candidates are weakly interacting massive particles (WIMPS) which were in thermal equilibrium with the species in the cosmic plasma via weak interactions [13]. In this scenario, galaxies form in a "bottom-up" fashion: small scales (favored by the small velocity dispersion of CDM particles) become non-linear and collapse first, and their merging and accretion leads to formation of structures on larger scales. On these scales, data of the structure of the universe is consistent with CDM driving the formation of galaxies and clusters, however, no viable fundamental particle within the standard model (SM) fulfills these properties [3,14] .
The standing ΛCDM paradigm, is in remarkable agreement with large scale cosmological observations (see for instance [15][16][17]) and it is also compatible with an increasing amount of observed galaxy properties (e.g. [18] and [19]). However, it has been noted that in this paradigm it is challenging to describe some observables on smaller scales, such as the "missing" dark matter sub-halos or the so called core-cusp discrepancy [20]. High resolution cosmological simulations of average-sized halos in ΛCDM predicts [21] an overproduction of small-scale structures, significantly larger that the observed number of small satellite galaxies in the Local Group [22,23]. Moreover, N-body simulations of CDM-only predict a singular density profile for virialized halos [24,25], while observational evidence points to dwarf spheroidal galaxies (dSphs) having smooth cores in their central regions [26,27]. Some other tensions have been raised between CDM-only predictions and observations (see for example a review in [20]).
Among the earliest approaches to alleviate/resolve those conflicts is to consider two DM components, one "cold" and one "hot" (C+HDM) [28,29]. More recent models feature only warm dark matter particles (WDM) [30], meaning that they are semi-relativistic during the earliest stages of structure formation with non-negligible free-streaming particle length. WDM models feature an intermediate velocity dispersion between HDM and CDM that results in a suppression of structures at small scales due to free-streaming [31]. If this free streaming scale today is smaller than the size of galaxy clusters, it can provide a solution to the missing satellites problem [26], as well as explain the smooth cores of dSph without spoiling the CDM predictions on large scales [32][33][34]. A particular, promising realization of these WDM models has been the minimal extension of the SM by intermediate-mass sterile neutrinos in the O(keV) range known as νMSM (see, for example, [35] for a review). From the astrophysical point of view, fermion masses in this range and up to O(0.1 MeV) seem to also be favored by recent elementary particle based DM halo studies [34,36], where selfgravitating equilibrium systems were shown to be both in excellent agreement with rotation curve observations while thermodynamically stable (coarse-grained entropy maxima) within cosmological timescales [37].
Another compelling alternative to colissionless CDM, apart from WDM, is to consider interactions in CDM. This consideration relaxes the assumption that CDM interacts only gravitationally after early decoupling, and includes interactions either between DM and SM particles or additional hidden particles, or among DM particles themselves. These later models are denominated as "self-interacting" DM models (SIDM) (see [38,39] for reviews). Born out of N-body simulations [40], SIDM halos could explain the cores of galaxies when a 2 ↔ 2 interaction is assumed, with cross-sections constrained to be roughly of σ/m ∼ 0.5 − 10 cm 2 /g [38,41,42]. However, certain tensions have been raised about the upper limits in the self-interaction cross section, motivated by a more refined analysis of the Bullet Cluster [43]. Likely indicating the need of a velocity dependent cross sections (i.e. σ as a function of the rms velocity of DM particles) sensitive to the baryonic environment [44].
Most SIDM studies assume a cosmological evolution identical to CDM, implying that the linear matter power spectrum remains unchanged. However, many models include other ingredients that can produce small scale damping by acoustic oscillations reminiscent of WDM [45][46][47][48]. A good example of the later are the DM + Dark Radiation (DR) models considered by the ETHOS collaboration [45], who created a framework for structure formation that encompasses several microphysical interaction models via an effective theory. Interestingly, interacting scenarios combining DM+DR interactions with SIDM effects, are able to generate a truncation in the power spectrum while producing shallower inner density profiles [49], alleviating the core-cusp and missing satellite problems altogether.
So far, we have mentioned both WDM and SIDM as possible solutions to the tensions between ΛCDM and observations on small scales, and discussed about their possible realizations. Here, we take both approaches into consideration. Previous studies have shown that the inclusion of self-interactions among WDM particles in quasi-relaxed DM halos can alleviate some constraints, as shown in [50,51] for the case of self-interacting right handed neutrinos. Also in [50] it is discussed the possibility of novel sterile neutrino production mechanisms through heavy mediators, while further effects of including a scalar self-interaction in the νMSM active-sterile mixing production scenarios, were considered in [52].
We focus here on the description and treatment of the linear theory of cosmological perturbations for self-interacting WDM (SI-WDM) scenarios, and provide explicit expressions for the Boltzmann hierarchies for different self-interacting sterile neutrino DM scenarios, with its corresponding beyond SM field mediators. Efforts on calculating the evolution of these perturbations either in traditional CDM or WDM scenarios (see [53] for a summary), have been outlined either via semi-analytic methods such as in [54,55], or via numerical integration of the coupled Einstein-Boltzmann system [53,[56][57][58]. For the latter, freely available numerical routines such as CAMB [57] or CLASS [58] exist as general purpose tools, or more specialized ones as the (CDM-based) ETHOS code [45] for interacting DM+DR models. An earlier work [59] pioneered the inclusion of SI-WDM on numerical Einstein-Boltzmann solvers (though under important simplifications, see also [60]), finding an enhanced suppression of power in small scales when compared to WDM only evolution.
The objective of this work is to contribute to the findings of these early realizations of SI-WDM structure formation. To this aim we provide here a systematic and accurate treatment of collisions in WDM models extending [59], and at the same time retaining certain independence from a particular Lagrangian self-interacting model. Our procedure is motivated by the tools provided by the Boltzmann hierarchies for interactive (active) neutrinos [61][62][63]. They are used and generalized to perform an accurate framework for the collision term in the linearized Boltzmann equation for the SI-WDM species, and derive an explicit and analytical expression for the equations of motion. Motivated by [45], we do not commit to a particular form of the scattering amplitude, but provide a general parametrization in terms of model dependent coefficients that naturally includes several interaction mediators such as a massive scalar (as seen in [52,62,63]) or a vector field (as proposed from first principles in [50,51]). The general results here presented are aimed (but not limited) to further evaluate the SI-WDM effects in the matter power spectrum, CMB anisotropies, halo models and production mechanisms, and may also be useful beyond the study of DM such as the study of active neutrino physics and their anomalies [63].
In order to set our notation and conventions, in what remains of this section we briefly introduce the cosmological perturbation theory and Einstein-Boltzmann equations.

Cosmological Perturbation Theory
In cosmology, the evolution of perturbations to the isotropic homogeneous background, which are originated through a primordial power spectrum and will eventually collapse to form the myriad of observed structures today, is handled through the Einstein equations. There, the universal spacetime metric is split into a background Friedmann-Robertson-Walker (FRW) metric and a small perturbation to said metric. The Einstein equations govern the evolution of this perturbation with the perturbed energy-momentum tensors acting as sources. Several choices exist in order to describe these metric perturbations: a "gauge freedom" in the equations. Here, we will use the so called synchronous gauge, where the line element is defined as where the scalar mode of the perturbation h ij can be described in terms of two fields h( k, τ ) and η( k, τ ) as A discussion on gauge freedom and gauge modes in the context of perturbations to the FRW metric can be found in [53,64,65]. Here, we quote the final form of the Einstein equations in the synchronous gauge, in Fourier space: with Σ the traceless component of T i j , ρ and P the background density and pressure respectively; and the metric perturbation functions h, η in Synchronous gauge are defined as in eq. (1.2).

The Relativistic Boltzmann Equation
In order to close the system of equations in (1.3) without the assumption of a perfect fluid, the perturbations in a given energy component can be obtained in a more general way by making use of the Boltzmann equation, which governs the evolution of the phase space distribution function (DF). As a relativistic invariant, this function is used to describe the number of particles of a given fluid in a differential unit of volume: where x and P are the spatial coordinates and their conjugate momentum respectively. In practice, it is more convenient to describe the perturbations to this function as a function of comoving proper momentum q i ≡ a(τ )p i (with p i measured in a comoving frame) as: where q = qn is the comoving momentum, and n j its j direction component and f 0 is the background DF. The phase space density evolves according to the relativistic Boltzmann equation. In terms of these new variables, this is: where Γ γ αβ is the general relativistic metric connection and P α refers to the components of the conjugate momentum, defined as P α ≡ (E/a, a[δ ij + h ij /2]p j ) in terms of the 4-momentum p j measured by an observer comoving with the FLRW coordinates and h ij the metric perturbation in the synchronous gauge (see [53] for details) . The right hand side of the equation involves the terms due to collisions (referred here as the collision term), whose form depends on the type of particle interactions involved. In the case of a general-relativistic formulation of perturbations, the derivatives with respect to the coordinates df /dx i and df /dq depend explicitly on the way one chooses to express the perturbed metric: the so called "gauge choice". In this case, we will use the synchronous gauge, where the metric takes the form of eq. (1.1).
We refer the reader to [64] for a comprehensive explanation on perturbed FRW metrics and the different gauge choices, and [53] for a "canonical" application to most of the cosmological fluids in more that one gauge. In k-space, the equation that dictates the evolution of the perturbation to the phase space distribution Ψ ≡ F/f 0 can be obtained from (1.7) and (1.6), to first order in F as: with = aE the comoving energy and h, η the potential functions describing the scalar mode of h ij defined as in eq. (1.2). This equation is to be solved together with the Einstein equation to give a closed system. These equations provide the evolution of the metric perturbations with the perturbations of the total energy-momentum tensor (built as the sum of the contributions of all relevant components) acting as a source term.

The Boltzmann Equation for SI-WDM: Interaction Terms
Here, we set our focus on the right-hand-side (RHS) of equation (1.8). This term describes the interaction between the different particle species, and the eventual self-interactions between the same species. Some components can be considered as collisionless during most of their lifetime such as CDM [53]: for some approaches on eventual collisions see e.g. [45,47,66,67]. For most other components the collision term plays a major role in their evolution, such as in the case for photons and baryons.
Several advancements have been made in the past years in the direction of dealing with the collision term in cosmological simulations from first principles. See for example [45,47] for a streamlining on the treatment of the term in CDM models, and [66,67] for an approximation of the full term by a Fokker-Planck operator. The focus of the following section is to extend the works of Oldengott et. al. [61,62], where the full term has been calculated for ultrarelativistic species (active neutrinos) and scalar mediators for the interactions between particles. Such an extension is centered in the case of SI-WDM, including a more general scattering amplitude for species that are neither ultrarelativistic nor fully non-relativistic at decoupling, with emphasis in self-interacting sterile neutrinos.
The RHS of equation (1.8) counts the number of collisions a particle species i undergoes in a time interval dτ per unit phase space. For a CP invariant two body scattering process i + j ↔ m + n, the full expression for the collision term is: where g is the number of internal degrees of freedom of each species, |M| 2 is the squared Feynmann amplitude for the process and δ (4) D is the Dirac delta functional over the energymomentum 4-vectors labeled with boldface.
The zero-order integral, which dictates the evolution of the background phase space distribution f 0 is simplified under the same assumptions as: (2. 2) The first order collision integral, which involves the first order perturbation F ( k, q, τ ) can be simplified in the case of interactions ii ↔ ii to: where we have made use of the symmetry of |M| 2 under the exchange q ↔ l , and under the assumption that Bose enhancement and Pauli blocking are negligible as is customary done for DM candidates on such early epochs [68]. In the case of a cosmological component that only interacts with itself, this would provide a source term in the RHS of equation (1.8), the equation of motion for the phase space perturbation F . The zero and first order Boltzmann equations form a set together that must be solved for each particle species under consideration, along with the Einstein equations for the evolution of metric perturbations.
It is typical on studies of self-interactions to assume a certain form for the background distribution functions [45,59,61] such as, for example, an equilibrium distribution which implicitly assumes either a thermal decoupling history of DM or a period of strong coupling in self-interactions. 1 Here, we will maintain a general form for f 0 and provide the necessary collision term to obtain its evolution via the zero-order Boltzmann equation.
In [61] both the first (2.1) and zero order collision integrals have been considered for the case of active neutrinos with a scalar interaction. In that case, a specific interaction model has been evaluated and the particle mass of the neutrinos has been neglected, given that they remain ultrarelativistic until late times. Here, we maintain a certain level of generality in the choice of interaction amplitudes, and explicitly include the particle's mass. This generalization of the collision term can be useful in certain WDM models that include self-interactions between dark particles. In particular, for those models, where the ultra relativistic to non relativistic transition takes place in the radiation dominated era, the corresponding limiting cases do not provide accurate results [69]. These topics are more thoroughly discussed in the following sections.
Besides the above mentioned assumptions for the collision terms in (2.2) and (2.3), we focus here on the case where the only relevant source of interaction the self-interaction among the Dark Matter particles themselves. However, if other interactions are relevant our results can be generalized by adding the corresponding collision term to the RHS of (1.8). Moreover, the evolution of the mediator fields should in principle be studied self-consistently. Nevertheless, in certain situations one can neglect the backreaction of those fields. For instance, in the case of very massive mediator particles this assumption is justified as the population of the mediators should be Boltzmann-suppressed at the times of interest. This is generically not true, however, in the case of a massless mediator: the contribution of the mediator population to the energy-momentum tensor may not be negligible and the interactions between these two components should be properly accounted for. Here we do not address the dynamics of mediator fields and restrict our analysis to the effects of the self-interactions of WDM. 1 A typical example for beyond SM neutrinos is (assuming relativistic decoupling of Self-Interactions), to set f0 ∝ exp(−q/T dec,0 ), where q is the comoving momentum and T dec,0 is the SI decoupling temperature today.
In [61] an extra normalization factor is included to provide a correction accounting for the effects of Fermi statistics in the number density.

Scattering Amplitude
Further assumptions enter the expression we will use for the spin-averaged scattering amplitude |M|. We will assume that this amplitude can be expressed as a second degree polynomial in (s, t): This parametrization leaves the m ij coefficients free as model dependent constants and allows us to recover a few relevant cases for our study, such as the ones to be considered in 3.3. This assumption, together with the approaches taken in describing both the collision terms and Boltzmann equations, allow us to complement previous works [45] aiming to describe selfinteracting species in cosmology. It is in this way that we maintain some model independence, being able to describe a wide array of (elastic) interaction cross sections either in an exact or approximate way.
Particularly, this parametrization recovers both the tree level self-interaction due to a massless mediator (in the ultra-relativistic limit, which turns out to be a constant scattering amplitude) and most tree level interactions due to massive mediators with m med m, where m is the DM mass and m med is the one of the mediator. Notably, this includes both of the examples studied in [61] as well as many more.
In order to explicitly perform the collision term integrations, it is necessary to recast this expression into their respective powers of t, which reads (with B t and C t trivial functions of s), as relevant in the case of the C 3 , C 2 integrals as demonstrated in appendix A. A similar expression works in the case of the C 1 integrals, this time involving u (having used the relation s 2 + t 2 + u 2 = 4m 2 , with A u , B u and C u simple functions of t as shown in appendix A):

The First Order Collision Integral
In this section we write down the final results of the first order collision term in (2.3), and refer to the reader to appendix A for the detailed derivation. In what follows we set g i = 2 in order to simplify the notation. If g i = 2 it is immediate to generalize our results. In terms of integrations in energies and Mandelstam variables, C 3 can be expressed as: with χ(s) defined as and {A t , B t , C t } given in eq. (2.5).
Here and in what follows we use the convention that all integrals run over the full range of the respective variables unless it is explicitly specified.
The calculations for the C 2 term are identical to the ones developed in A.1 for C 3 . The only difference is that the roles of the background and perturbed DF are reversed. This can easily be seen from the definition of the term in (2.3). So, the final expression for the integral is where we have implicitly used that l is a function of only (E l , s). Given q, this is straightforward to check from the definitions of s, E l . In the case of C 1 , the calculation diverges greatly from the one of C 3 . In this case, both the background DF and the perturbation are integrated over, and to perform the integration it is necessary to know f 0 (E q , τ ). This integral can be expressed in terms of time-dependent collision kernel K(E q , E q , t, τ ) as put in terms of moments of the background distribution function f 0 , that take the form which are functions of (E q , E q , t) only through R 2 , defined as This kernel is the most complex part of the collision term, mainly due to its explicit dependence on time through the momenta of the background DF. However, once the scattering amplitude |M| 2 is specified, it should be numerically feasible to evaluate the integrals.

The Zero Order Collision Integral
The treatment of the term D 2 [f ] mimics exactly the one for C 3 [f ] but with the simplification F ( q, k, τ ) → f 0 (E q ). Thus, this term can be expressed as The term D 1 [f ] is much more complicated. The key to solve this integral is to define a method to recast an integration in an angular variable by an integral in energy, as described in appendix B.1. As shown in such appendix, this procedure leads to a collision integral that can be expressed as in terms of 4 integrals of kernel functions: the integration limits I i are defined in (B.7) schematically as and the kernels k i (E q , E q , E l , τ ) are given by 9) and t 1,2 defined as the two solutions to the following equation: Thus, we have arrived at a somewhat general expression for the collision integrals which depends on the scattering amplitude only through the coefficients A t,u , B t,u and C t,u defined in equations (2.5) and (2.6). After this, in order to obtain a Boltzmann hierarchy that can be in principle solved numerically one can follow a procedure analogous to the one described in [61]. We perform this procedure in section 4, while next we provide some examples of the kernel functions obtained from different models of the self-interaction.

Kernel Functions for Different Mediator Models
In this section, we calculate the different kernel functions involved in the collision integrals for a small subset of self-interaction models. We need to compute the coefficients Motivated by the possibility that the DM constituents are sterile neutrinos, we consider the following three cases: the first two evaluated by [61] which are interactions mediated by scalar particles in both the very massive and strictly massless limits; and the the case of a heavy vector field proposed in [70].
For the first two cases, the interaction Lagrangian can be written as (further information about the scattering processes can be found in [61]): where g is the scalar coupling constant, φ is the scalar field and ν R is the DM field modelled as a right handed neutrino. In the case studied in [61] the massless scalar limit reduces to a constant amplitude, however, this does not happen generally. We refer to [71] for an expression of the scattering amplitude for scalar mediators of arbitrary mass. In this study, as an example, we only consider a constant amplitude case, reminding that only in the massless limit it corresponds to a zero mass scalar mediator. In the second model we consider a massive scalar mediator, assuming that m med = m φ E (with E denoting the mean energy of the colliding DM particles).
The vectorial model of [70] also assumes DM is given by right handed neutrinos but with an interaction Lagrangian given by with g V acting as a coupling constant and V µ , the massive vector field. In the cases considered here, and under the assumption g V 1 all mediators fall into the massive case m med = m V E (see (3.25)).
The authors of ref. [70] have proposed this interaction to derive an effective selfinteraction for a self-gravitating system of sterile neutrinos on galaxy scales. It was there shown how a O(10 1 ) keV-fermionic DM concentration at the center of the DM halo of the Milky Way, could work as an alternative to the massive black hole (SMBH) in SgrA*. At the same time the model may work as an alternative DM candidate that could provide a different explanation to the small scale cosmological observables.
It is important to note that, while motivated by the study of sterile neutrinos, these expressions and interaction models remain general and can be used in other applications such as (massive) active neutrino cosmology.

Constant Amplitude
We start with the simplest of models: a constant scattering amplitude |M| 2 . We adopt the notation used in [61] for the massless scalar mediator. This constant amplitude can be expressed as where g is the scalar coupling constant in the ultra-relativistic case. Being constant in the involved momenta, the coefficients of the expansion in Mandelstam variables are, simply: For the χ(s) function appearing in the final form for C 3 and C 2 , we obtain For the time dependent kernel function K(E q , E q , t, τ ) in C 1 , we find with f 0 0 as defined in eq. (3.5) and with t 1,2 defined as the roots of the equation P (t) = 2E l − E q + E q with t 2 > t 1 , as used in section 3.2 for the kernel functions to calculate the background DF in D 1 . When comparing these expressions with the ones used in [61], both collision integrals coincide in the limit m → 0, showing explicitly that our more general expression for the collision term reduces to this known limiting case.

Massive Scalar Mediator
We follow here the considerations of [61] for the case of a scalar mediator which is considerable more massive than the scatterers. In this case the population of scalar particles would be Boltzmann suppressed, so there would be no need to track the evolution of their population. Moreover, in this scenario, the neutrinos would be initially in thermal equilibrium (as noted in [61] and references therein). In this case, the interaction amplitude reduces to: with m φ denoting the scalar mediator mass. Here, by using the identity s + t + u = 4m 2 we can either replace u or s in the scattering amplitude to obtain the two sets of scattering coefficients [A, B, C]: .
Making use of these coefficients, the kernel functions χ, k i and K as defined in section 3.2 and 3.1 respectively, read as follows: 22) where the integration regions in the t variable for these last integrals are in the range Here the background DF kernels k i have analytical (though complicated) expressions in (E q , E q , E l , t 1,2 ), which are not very illuminating to write down.

Massive Vector Field
In [70] the authors calculate the 4-fermion self-scattering amplitude for the right handed sterile neutrinos with the interaction (3.11), and reach the following result: where m V is the mediator mass, θ W the (dark sector) Weinberg angle and pq ≡ p µ q µ on a 4-vector notation. By making use of the following Mandelstam variables properties the equation (3.23) can be put in terms of (s, t, u). As for massive scalars, the scattering coefficients can be calculated using s + t + u = 4m 2 : , (3.26) and the expressions for the kernels are: 29) where, again, the integration regions in the t variable for these last integrals are in the range [t cos θ=−1 , t cos θ=1 ] for i = 2, 4 and [t 1 , t 2 ], where t 1,2 are the solutions of P (t) = 2E l − E q + E q for i = 1, 3. These last kernel functions for the background DF have analytical forms but are not illuminating, just as in the massive scalar case (see appendix A for details).

Boltzmann Hierarchy
Once having obtained the expressions for the collision integral kernels, it is a standard practice to perform a Legendre expansion in (1.8) in order to construct a so called Boltzmann hierarchy of equations, which is independent of the angle between k and q. In order to calculate the time dependent kernels above, the full solution to the background DF f 0 (E q , τ ) must be obtained. If we assume that, in the time scales of interest, the collision term is only due to self-scattering, then the evolution of f 0 is governed by: Here, we will assume that the initial conditions for f 0 are set beforehand at some early time and that its subsequent evolution is only governed by self-interactions as said before. This ansatz, notably, excludes the production mechanism that should give rise to the initial population of these particles: we follow [59] and implicitly assume that the mechanism for production is not significanly affected by the self-interaction mechanism. Once the evolution of this distribution is known, the various moments f 0 i defined in eq. (A.68) can be calculated, thus allowing to obtain the time dependent kernels K.
As the coefficients in the LHS of the Boltzmann eq. (1.8) are only functions of | q|, | k| and cos ≡k.q, it is assumed that the RHS also depends only on these parameters. Thus, in order to express the relevant equations of both sides in a Legendre series we first write where P l (cos ) is the l-th Legendre polynomial and F l is the l-th multiple of the perturbed DF.
In this case however, a residual dependence on the azimuthal angle between q and k, ψ, is still present. In [61] the authors argue that an averaging over ψ in the collision terms has no effect on the expressions for the collision integrals, and they perform such average. The argument is based on the fact that the LHS of (1.8) is not affected by such averaging. Also, from a phenomenological viewpoint, the only observable is the integrated effect of the perturbation, further strengthening the claim. Therefore, in what follows, we perform such average as well.
The moment decomposition in the LHS of the Boltzmann equation is well known, and the reader can consult [67] for the expressions for massive neutrinos in the collisionless case on both typical gauge choices.
In order to calculate the moment expansions of the collision integral, we make use of the following property: given a collision term with the form ∂f ∂τ the (ψ averaged) l-th multipole can be written as: However, our expressions for the collision integrals are expressed in terms of Mandelstam variables and energies instead of angles and momenta. Our kernels are also expressed in terms of these variables. We can solve both problems by recasting the integration in Mandelstam variables in the definition of the kernel moments. So, for a collision integral of the form ∂f ∂τ the property (4.6) would be modified as follows: with K l (E q , E l , τ ) ≡ dsK(E q , E l , s, τ )P l (cos θ(s)) . (4.10) Then, putting together the results of sections A.1, A.2 and A.3 we arrive at the following moment expansion, in Synchronous gauge: with the various kernel moments defined as: where K l , χ l are the Legendre transforms of the K, χ kernel functions defined as in (4.10), and where the l-th moment of the perturbed DF F l is defined as in (4.4), and we have chosen to express the momentum dependence in terms of energy for consistency. In order to solve this hierarchy, the kernel functions for the interaction model must be specified.

Relaxation Time Approximation
Even if the evolution of the background DF f 0 may seem complicated due to the collisions in play, its effect may be accounted for in a much simpler way depending on the particularities of the interaction. Concretely, if the rate of particle interactions is much higher than the rate of expansion of the universe, measured roughly by H, the Hubble rate, we may assume that the shape of the distribution function is one in equilibrium. That is to say, a DF that obeys such as Maxwell-Boltzmann, Fermi-Dirac or Bose-Einstein distributions, depending on the particle model used. It is possible to construct a substitute collision operator for f 0 that reconstructs the expected behavior for small departures from thermal equilibrium: This is know as the relaxation time approximation of the collision operator. The relaxation time τ is the timescale in which the system is expected to relax to equilibrium. This parameter can in principle have a q dependence and is commonly defined as [72]: which involves the integral of the kernel χ, evaluated in the thermal equilibrium background DF. It is straighforward to evaluate these integrals, as they only involve known functions. We have performed these numerically for the interaction models posed in 3.3, and the results can be seen in figure 1 (we refer to appendix C for more details).
For all of these models, we follow [59] and assume that the abundance of WDM and its primordial distribution function are already set deep into the radiation dominated epoch and the effects of self-interactions in these initial conditions can be effectively decoupled from the evolution of perturbations.

Application to Self-Interaction Decoupling
In order to evaluate whether or not a thermal f 0 can be assumed, for a given particle physics model for the interaction, one may look at the ensemble averages of the interaction rate Γ(E q ) = τ −1 (E q ). This value is to be compared to H at this point: if Γ H, the system is effectively in thermal equilibrium and adopts an equilibrium background distribution function f 0 .  [70]). The calculations are performed for a DM particle mass of 10 keV, and the vertical line marks the relativistic-nonrelativistic transition temperature. This is equivalent to the traditional approach used to determine if a species has decoupled from the rest of the cosmic plasma in the standard sector (see for example [62]). In other words, a given interaction is considered to cease being relevant if the interaction rate per particle Γ ∼ τ −1 th (see appendix C) is overtaken by the Hubble expansion rate, which in the radiation dominated era is H ∼ T 2 /m pl where m pl is the Planck mass. So, in summary, a thermal background DF can be assumed if at some point during the evolution of the perturbations, the self-interactions were a dominant phenomenon in the sense Γ > H.
We can see how this interaction rate evolves along with the temperature of the plasma in fig. 2, using the models of section 3.3. We work here on the assumption that the self-interaction decouples (this is to say, Γ > H) while the particle itself is still relativistic. At the moment the self-interaction decouples its distribution function remains "frozen-out": the function itself remains unchanged and the evolution is just given by the redshift in physical momenta p ∝ a. If decoupled while relativistic but well after the initial production of these particles, the distribution is frozen with a form f 0 ∝ e −p dec /T dec and the redshift in momenta can be reinterpreted as a temperature evolution of the form T ∝ a −1 . In the case of a non-relativistic decoupling of the self-interaction however, the distribution is frozen out at f 0 ∝ e −p 2 dec /(2mT dec ) and the temperature is interpreted to evolve as T ∝ a −2 .
We can see in fig. 2 that the assumption of relativistic decoupling of the self-interaction does not necessarily hold for some interaction constants (e.g. C V ∼ 10 8 G f ). In that case, the self-interactions should decouple while non-relativistic, if at all, and may alter the distribution function which in turn may render the method we used to obtain fig. 2 (described in (C.11)) inapplicable, as it assumes f 0 ∝ e −E/T . In particular, couplings C V ∼ 10 8 G f in the vector field case were shown in [70] to correspond to cross sections in the range σ/m ∼ 0.1−1 cm 2 /g, which are usually considered to alleviate various problems in N-body simulations on self-interacting CDM, and are strongly constrained by observations [38]. 2

An Approximate Form for the Collision Integrals
The full form of the Boltzmann hierarchy for these species (4.11) can be reduced by making use of the relaxation time approximation. The most straightforward way to do this is by simply replacing expression (5.2) into the collision term and calculating the new hierarchies, through expression D 2 [f 0 ] (3.6). In [59] such an approach is taken in a simplified way: instead of the full (energy dependent) relaxation time τ rel , its thermal average is used (see appendix C): This simple approach however leads to an important conceptual error: this approximation (and to a certain extent, also form (5.2)) qualitatively simply "erase" the perturbations to the DF F ( q, k, τ ) [73]. This violates conservation of particle, momentum and energy densities, resulting in a poor approximation to the full collision term in the case of perturbations. In [59] it is noted that this can be avoided by setting the C l=0,1 [f ] = 0, and these conservation laws are recovered, thus arriving to a Boltzmann hierarchy of the form: This relaxation time approximation [59,62,74] has the advantage of localizing the equations in momenta, which results in more efficient numerical integration by eliminating all coupling between different momentum bins and allowing for sparse evaluation. This simplified approach can be generalized slightly by using instead the separable ansatz namely, assuming that the "temperature perturbation" F l (E q , k, τ ) is independent of momentum, that is Then, the l-th collision term can be reduced to: with: While this is a significant simplification to the collision term, as proposed in [62], further simplifications can be done by performing a momentum average Under this approximation the Boltzmann hierarchy reduces to 3 : 10) provided C l=0,1 = 0. This approximation can be further reduced to the form (5.5) by assuming α l≥2 = 1. If it results on a better overall approximation than (5.5), remains to be explored in future works.

Summary and Outlook
Throughout this work, we aimed to fill a gap in the description and treatment of linear theory of scalar perturbations in cosmology by including the case of a self-interacting warm dark component. Motivated by the possible impact of these self-interactions in large and small structure formation scales, we provide an accurate treatment of collisions in the early universe, extending previous works on the subject while maintaining a phenomenological approach that allows us to retain certain model independence on the particular interaction Lagrangian. By extending the treatment in [61][62][63] for active neutrinos, we calculated the first and zero order collision terms and provided a general framework in order to include these collision terms in the coupled Einstein-Boltzmann system. This was done with the objective of accurately evaluating the effect of WDM self-interactions on the linear power spectrum and the CMB anisotropies.
In section 2 we provide a short summary about the assumptions used in this framework, as well as the general form (2.4) for the interaction amplitude that was used. It is shown there that this form can exactly describe several models of massive mediator interactions between sterile neutrinos, though not limited to those cases, and including for Majoron-like scalar mediators between right handed neutrinos. The main calculations are given in section 3.1 and 3.2 where we show the results for the first and zero order collision terms respectively in the SI-WDM scenario. Also, we provide in section 3.3 some examples for a handful of specific interaction models, along with the corresponding coefficients for the collision terms. A detailed treatment on how to include these collision terms in a Boltzmann hierarchy is shown in section 4, and several possible simplifications in order to treat the evolution of the background and perturbed distributions based on the relaxation time approximation are shown in section 5.
While developed with the intent of being used in the calculation of cosmological perturbations in the case of SI-WDM, the forms of the collision terms themselves are quite general and they can be used also in several other applications, for example in the case of (massive) active neutrino cosmology, as discussed in section 3.3. The implementation of this formalism in a CMB Boltzmann solver in order to quantitatively explore the effects of the SI-WDM framework in cosmology, is left for future work (in progress). In this direction, further exploration of possible approximation schemes or efficient computation methods is key in order to successfully implement the hierarchies (4.11) in a practical way. Moreover, the phenomenological assumption of (2.4) for the interaction amplitude can be expanded upon: a more general form can possibly include other types of interaction Lagrangians, such as light mediators or more complex models accurately. Also, it should be possible to extend this formalism to include the effects of Bose enhancement or Pauli blocking by generalizing the collision kernels, as was argued in [61].

A.1 Calculation of C 3 [f ]
The term C 3 can be expressed as: where g i is the number of degrees of freedom of the particle and we have omitted the k, τ dependencies on the distribution functions for compactness.

A.1.1 Solving for l and β
The energy conservation Dirac delta δ (4) D (q + l − q − l ) can be used to solve directly the l integral using the definition of the invariant integration measure: and Θ is the Heaviside theta function. The following parametrization is used for the momentum 3-vectors: So the argument of the Dirac delta in (A.3) can be expressed as: In this parametrization for the momentum 3-vectors, the integrals in the collision term can be expressed as g( q, l, q )) , where we have used d 3 l = 2πl 2 dld(cos α) and d 3 q = dβq 2 dq d(cos θ), we have omitted the relevant integration bounds except on the β integral and assumed that the scattering amplitude |M| does not depend on the azimuthal angle β. For this integral, we use the following identity of the Dirac delta: that allows calculations of compositions between delta functionals and functions, where β i are real roots of the real function g(..., β, ...). For this integral, we have From this, we can infer that two solutions β i exist: one in the interval [0, π] and one in [π, 2π].
As the absolute value of the derivative of g is the same in both solutions, we can express the integral in (A.6) as To ensure the physical condition that | cos β i | < 1, we add a Heaviside step function in cos 2 β i . The following property follows from (A.9) : So, the C 3 integral in (A.6) can be expressed as Here, an important point in the calculation is reached. The remaining angular integrals are in cos θ and cos α, the angles between q and q , and q and l, respectively. Now, without any knowledge of the background function f 0 , up to three of the remaining four integrals could be solved. However to continue solving from here on we need to know about the scattering amplitude M. We will go as far as possible without specifying this, and then we will assume an ansatz for a general form of M. To continue we express the term |∂g/∂β| cos β i in terms of the variables {q, l, θ, α} as follows: with coefficients: (A.14) Now, let us consider in greater detail the integration of the Heaviside theta function in equation (A.12). The argument of the function is a quadratic function in cos θ with a negative leading coefficient. Thus, the function will only be non zero if two real roots of the polynomial |∂g/∂β| 2 cos β i exist and it will be unity between them. So, the Θ function can be translated into a border condition for the cos θ integral and an existence condition for the roots: with x 1,2 the roots of the polynomial |∂g/∂β| 2 cos β i , The argument of the Heaviside step function can be expressed as cos y) , So this argument has 4 real roots in E q : {−m, +m, R 1 , R 2 }, with R 1,2 given by The −m root can already be discarded as non physical. In order to to obtain the ordering for the rest of the roots, thus the non-zero intervals for the Heaviside step functions, we will next develop some alternative notation for the angular variables θ, α, making use of Mandelstam variables.

A.1.2 Mandelstam variables for C 3
Mandesltam variables are Lorentz invariant quantities constructed with the relevant information on a two on two scattering process and are defined as: We can make use of these quantities to advance in the remaining integrals for the collision integral C 3 . First of all it is important to note that at this stage the derivations here and in [61] start to diverge. Some of the interesting properties of Mandelstam variables, which allows them to be of use in these calculations, are related to the center of momentum (CoM) frame that cannot be properly defined in the case of collisions between massless particles. Indeed, for the case of identical particles the Mandelstam variables can be calculated in the CoM frame making use of their Lorentz invariance: with E CoM m , p CoM m the individual particle's energy and momentum magnitude measured in the CoM frame and θ CoM the scattering angle measured in the same frame. We can also express these variables in terms of the quantities we have used throughout the calculation of C 3 as: where here, the quantities are measured in the "lab" frame, that is to say, the fixed frame in which we have measured q. It would not be possible to change the whole integral to CoM quantities making use of the Lorentz invariant measure d 3 p/2E p because q is fixed by the LHS of the Bolztmann equation. As a first use for these quantities, let us consider the roots R 1,2 in the polynomial above. The roots can be immediately recast as a function of s: From expression (A.21) we can see that s ≥ 4m 2 , thus these roots are always real. As for the ordering of the roots, it should be quite obvious that R 2 > m, R 1 . It is also possible to prove that R 1 > m: • Let us first define the total momentum 4 vector: p µ = (E q + E l , q + l). If we recall the form of a Lorentz transformation over a 4-vector, we can infer that the form of R 1 is very similar to the 0 component of this 4-vector, measured in a different coordinate system (Lorentz boosted). The magnitude of this boost from the "laboratory" system ( where we measure q, l ) to this new system would be | β| = 1 − 4m 2 /s in the direction of the total 3-momentum q + l 4 . The total energy p 0 measured in this new frame of reference is • On the other hand this energy, measured in the boosted system, can be also written as where, this time, m p µ = |p µ p µ | and γ is the Lorentz factor corresponding to the boost from the center of momentum system to this new system of reference.
• Equating these two expressions, one finds: 4 This is a valid boost rapidity β, as using the properties of the Mandelstam variables one can find , with these quantities being the individual particle's energies and momenta measured in the center of momentum system. So, going back to the expression for the last Heaviside theta (A.17) we can conclude that it is both non zero and physical only between R 1 and R 2 . We can then express the whole C 3 integral as where q (R) = √ R 2 − m 2 and we have made use of the fact that R 2 ≥ E q + E l to eliminate Θ(E q + E l − E q ). The observant reader may have noticed that not only we omitted the dependency on the dynamical variables of |M| 2 , but we have also simply carried it outside of both integrations, on l and β. As the Mandelstam variables encompass all of the relevant invariant quantities involved on the kinematics of the process itself, it is reasonable to expect the scattering amplitude |M| 2 to only depend on (s, t, u). Now, we have solved first a three momentum integral in l imposing the momentum conservation Dirac delta. When calculating scattering amplitudes momentum conservation is explicitly imposed, so we can consider |M| to already be evaluated at l = − q − l + q . Now for the β integral we have solved a Dirac delta in the function g defined in (A.3). If we express this function in Mandelstam variables we can see that it is simply δ D (s + t + u − 4m 2 ). Thus, we can impose the condition s + t + u = 4m 2 in |M| and bring it out of the β integral. Note that this condition is trivially fulfilled given the definition of s, t and u.

A.1.3 Change of variables to {s, t, E q , E q }
At this point, as mentioned above, further integration is not possible without specific knowledge of the scattering amplitude |M|. However in terms of the variables we have used so far, namely { q, q , l, l }, the scattering amplitude may indeed have a very complicated and ultimately redundant expression. The Mandelstam variables contain all of the Lorentz invariant quantities that are involved in the expression of |M| so that its only dependencies would be on (s, t), after having applied the identity s + t + u = 4m 2 . In order to accommodate to a more general expression for the scattering amplitude, it is convenient to change variables in the integrals of C 3 to obtain an expression involving (s, t). Indeed, we can make the following variable change: (A.29) The integration measures change accordingly as: So that the whole integral can be expressed as Now, the expression (A.13) that states |∂g/∂β| 2 cos β i as a second degree polynomial in cos θ can be recast into a polynomial in t using (A.22): The ingredients are all set to perform the t integral given an expression for |M|. We then assume that this amplitude can be expressed as a second degree polynomial in (s, t), as in (2.4): In order to explicitly perform the t integral, we group the coefficients in (A.34) in their respective powers of t: So the t integral becomes .
Then, we use the following identity: where x (+,−) are the solutions of ax 2 + bx + c = 0 to arrive at The expressions in brackets have fairly complicated forms, but mostly polynomial in nature in (E q , E q , E l , s, m). They can be factored into powers of E q in order to integrate them as: Once these are replaced into (A.37), we can perform the integral in E q using In this way we arrive at the following expression for the integral C 3 : which can be recast into a more compact notation: where the function χ is defined in (3.2), and we explicitly wrote the arguments of each function. Now, in this particular integral we could in fact formally integrate in s to obtain only a last integral over E l of the background DF times an arbitrary function. We leave this as it is, because this form will become practical in calculating the next collision integral C 2 . It becomes quite remarkable that the χ function only depends on s (which in turn is only proportional to the CoM energy). This might suggest that a more straightforward method of obtaining these integrals may be available.

A.2 Calculation of C 2 [f ]
The calculations for the C 2 term are identical to the ones developed in the previous section. The only difference here is that the roles of the background and perturbed DF are reversed. This can easily be seen from the definition of the term in (2.3). So, the final expression for the integral is where we have implicitly used l ≡ l(E l , s). This is true given the parametrization (A.4) for this vector, as its magnitude can be uniquely determined by E l and its angle with the vector q by the Mandelstam variable s. Therefore, unlike the previous case, the integral on s cannot be performed, since F also depends on it.

A.3 Calculation of C 1 [f ]
The procedure for calculating C 1 are very similar to the ones used to calculate C 3 , but with some changes in the parametrization of the vectors involved. Following from the expression (2.3), the full integral can be expressed as

A.3.1 Solving for l and β
Since the integrand does not depend on l, we perform the first integral using the momentum conservation Dirac delta along this variable (instead of l as it was used before). After doing this, we use a new parametrization that reflects our choice of remaining variables: Following the steps on section A.1, we reach the following expression: where this time the argument of the remaining Dirac delta g can be expressed as  We also add here a Heaviside theta function in cos 2 β i to ensure the condition | cos β i | < 1, but this time with the argument ∂g ∂β with coefficients: (A.52) Then the C 1 integral results: The argument for the last theta function is a second degree polynomial, this time in cos α, with negative leading coefficient. As before, we solve this condition by imposing integration limits in the cos α integral and adding a new Heaviside theta function to ensure the existence of real roots for the polynomial (A.51). The discriminant of this polynomial becomes: The argument of the function again has four real roots in E l : {−m, +m, R 1 , R 2 } with R 1,2 given by As before, we can discard any roots smaller than m as non physical. This time however, the ordering of these roots is different, as we will see below once we express these roots in terms of Mandelstam variables.

A.3.2 Mandelstam variables for C 1
So again, in order to assess the ordering of the roots R 1,2 , as well as facilitating the integration of the collision kernels, we make use of the Mandelstam variables. The relevant variables for the integrals here are t and u. Their expressions are given in (A.20), however as we have changed the parametrization of the vectors themselves for this collision term, we need their expressions in terms of the momentum variables: (A.57) Now, the roots of b 1 in terms of these new variables is We can see that the argument in the square root is always positive, because t is always negative. To find out the order of the roots, we start by noticing the following three conditions: 1. Thanks to the fact that t is spacelike, we can infer that |E q − E q | < | q − q |.
3. Thus R 1 ≤ 0 for any ordering of E q , E q .
Then, as R 1 is negative, it is not a physical value for the energy integral. The proof that R 2 ≥ m 2 , hence R 2 is physical, is slightly more complicated but possible, and is left to the reader. It is important to note in this case the ordering of these roots in order to understand the allowed ranges for E l . If we focus on the discriminant (A.54) we see that, as before, it has four roots but in this case the value is positive in the limit E l → ∞. So, as the ordering of all four simple roots is in this case [(−m ≤ R 1 ) or (R 1 ≤ −m)] ≤ m ≤ R 2 , we can conclude that the only range of energies where both E l ≥ m and the discriminant is positive is E l ≥ R 2 .
So, for C 1 we have boundary conditions for E l that are completely different from the ones for C 2 and C 3 . We can then use these limits in energy for the full integral but first we need to change variables to Mandelstam variables to perform them, these time using (t, u).

A.3.3 Change of variables to
So, we now make the appropriate change of variables proceeding in a very similar way as in the case of C 3 . In fact, the expression for the variables and differentials in E q and t are exactly the same as (A.29), (A.30). It is therefore immediate to see that the change of variables to E l and u is expressed as: with the corresponding integration measures given by Now, we had expressed |M| 2 in terms of only (s, t) in (A.34). However, this time s is not a relevant variable of integration. We can use the relation s + t + u = 4m 2 to obtain an analogous expression in variables (t, u): with the {A u , B u , C u } coefficients only depending on t. So, the C 1 collision integral can be expressed as: with y 1,2 the roots of |∂g/∂β| 2 cos β i in cos α. This expression depends on |∂g/∂β| 2 cos β i , this time given by (A.51), which can be recast into a polynomial in u as: with coefficients: (A.64) The integral in u can then be computed using the property (A.41), so that the full collision integral becomes The expressions for the coefficients accompanying {A u , B u } can be expressed in terms of powers of E l , in order to facilitate integration: (A.67) Here the approach diverges greatly from the one we took for C 2 and C 3 . In the previous cases we managed to also perform the integration in energy, but here this cannot be done without knowing the background DF f 0 . Since given f 0 , the integration in E l can be performed, it is convenient to define which is a function of (E q , E q , t) only through R 2 , and express the full collision integral in terms of these moments of f 0 : We can recast this expression in terms of a (fairly complex but mostly polynomial) integration kernel K(E q , E q , t, τ ).
This integration kernel is the most complex part of the collision term, mainly due to the explicit dependence of the integration kernel on time through the momenta of the background DF.

B The Zero Order Collision Term
In this appendix we provide details on the calculation of the zero-order collision integral. As shown in (2.2) this can be split into two parts: D 1 [f ] and D 2 [f ]. The treatment of the term D 2 [f ] mimics exactly the one for C 3 [f ] but with the simplification F ( q, k, τ ) → f 0 (E q ). Thus, it is immediate to see that this term can be expressed as in (3.6).

B.1 Calculation of D 1 [f ]
The term D 1 [f ] holds some similarity to the term C 1 , as it involves the integration of both of the vectors q , l on which the distribution functions are evaluated. The derivation for this term closely follows the one for the perturbed DF until equation (A.65). In this case, the result can be written as Continuing the procedure as we did for C 1 will require the moments for the background distribution function, which is exactly the quantity we are trying to obtain. Here instead of doing that, we will express this integral in such a way that it depends on the integral of a function times the background DF itself: we want to obtain a convenient expression for the integration kernel to be used in the equation for f 0 . To do that, we will integrate in t, given that the background DF does not depend explicitly on this variable. A complication in this approach comes from the fact that a simple modification of the order of integrals will not work, as the boundary for the integral in E l , R 2 , depends explicitly on t. In order to circumvent this, let us first define the variable P We change variables from t to this new variable P in the integral: Now, we rewrite the integration limits as a series of Heaviside theta functions Eq+E q |Eq−E q |+2m dP ∞ R 2 dE l = dP dE l Θ(E l − R 2 )Θ(P − (|E q − E q | + 2m))Θ(E q + E q − P ) .

(B.4)
Then, using the relations Here, each of the four parts is first integrated in an angular variable and then in the energies. Using this procedure for the integral (B.1) the collision integral can be expressed as: where I i refers to the integration method in energies specified in (B.7) and the kernel functions k i are expressed as in (3.8).
There, we have labeled the kernels with i = 1, ..., 4 according to the integration regions I i in the order in which they appear in equation (B.7). We have chosen to return the integrals to their original angular variable t instead of P , as the integrands themselves are functions of t and can't be expressed neatly in the new variable and the integration measure is also a function of both P and t. The integrals labeled i = 2, 4 possess the original limits on t, as the integration scheme has not modified the bounds in P appearing in (B.3). The ones labeled i = 1, 3 are the (two) solutions of the equation P (t 1,2 ) = 2E l − E q + E q . The integrands themselves are expressed in equation (3.9) on a similar way as in equation (A.69).

C Numerical Integration of D 2 [f eq 0 ]
In this section, we will provide the numerical integration scheme used to calculate τ . As was justified in 5, the expression for this quantity can be set in terms of the kernels introduced in 3.3: As we are not considering Bose enhancement or Pauli blocking, we set the equilibrium DF as a relativistic Maxwell-Boltzmann distribution f eq 0 (E q , t) ≡ f M B (E q , T (t)) = e −Eq/T , (C.1) and the function χ(s) depends of the particular interaction model used. Now, in order to numerically integrate these equations, it is necessary to adimensionalize this expression. We define:τ dΩσv mol E q E l (2π) 4 2 = π 4 d cos Θ cm |M| 2 1 − 4m 2 s , (C. 10) with Θ cm the scattering angle measured in the CoM system. This property, apart from relating χ and σ, reproduces the feature that σv mol E q E l is a function of s only (see [75]). We can use the procedures outlined in the calculation of thermal averaged annihilation cross section to obtain an alternative expression for the thermal average of the relaxation time: where T is the temperature of the background Maxwell-Boltzmann DF, K n is the modified Bessel function of the second kind of order n and N is a normalization constant defined as This expression yields the same results as the one in (C.8), but it is much easier and faster to implement numerically.