Momentum distributions of cosmic relics: Improved analysis

We solve coupled momentum-dependent Boltzmann equations for the phase space distribution of cosmic relic particles, without resorting to approximations of assuming kinetic equilibrium or neglecting backscattering or elastic interactions. Our method is amendable to precision numerical computations. To test it, we consider two benchmark models where the momentum dependence of dark matter distribution function is potentially important: a real singlet scalar extension near the Higgs resonance and a sterile neutrino dark matter model with a singlet scalar mediator. The singlet scalar example shows that the kinetic equilibrium may hold surprisingly well even near sharp resonances. However, the integrated method may underestimate the relic density by up to 40% in extreme cases. In the sterile neutrino dark matter model, we studied how the inclusion of previously ignored elastic interactions and processes with initial state sterile neutrinos could affect the nonthermal nature of their resulting distributions. Here the effects turned out to be negligible, proving the robustness of the earlier predictions.


I. INTRODUCTION
Current cosmological observations can be accommodated within the cold dark matter (CDM) paradigm [1]. This hypothesis is appealing within our present understanding of the structure of ordinary matter: extending the Standard Model (SM) of elementary particles and their interactions with dark matter particle degrees of freedom allows the abundance of CDM to be created by thermal production and decoupling of the dark matter particles in the expanding early Universe. The standard treatment for calculating this abundance of CDM particles relies on the Zel'dovich-Okun-Pikelner-Lee-Weinberg (ZOPLW) equation [2,3].
The ZOPLW equation is obtained from the Boltzmann equation by integrating over the phase space of the dark matter under the assumptions of detailed balance and kinetic equilibrium distributions [4][5][6]. However, the momentum distribution of dark matter may contain essential information that is neglected in this treatment. For example, if the dark matter production takes place at a resonance region, where the DM annihilation rate is strongly momentum dependent, the elastic reactions might not be fast enough to keep kinetic equilibrium. In this setting the true annihilation rate and hence the final DM abundance may deviate from the value obtained under the equilibrium assumption [6][7][8]. Another example concerns warm dark matter (WDM) whose momentum distribution may directly influence the cosmic structure formation by reducing the number of DM halos at small scales compared to CDM. In kinetic equilibrium, the suppression of the matter power spectrum can be well approximated via a single scale given by the WDM mass [9]. * kalle.ala-mattinen@helsinki.fi † matti.heikinheimo@helsinki.fi ‡ kimmo.kainulainen@jyu.fi § kimmo.i.tuominen@helsinki.fi However, if the DM particle is not in kinetic equilibrium, the resulting suppression may be more complicated [9][10][11][12]. Earlier calculations accounting for the DM momentum distributions tend to rely on simplifying approximations. For example, in the analysis of [10,11], the elastic interactions of the initial state DM particles have been neglected. In this paper, we complement these earlier analyses by presenting a numerical method based on discretization of the momentum space that allows for a completely general solution of the Boltzmann equations for the momentum distributions of multiple number of particle species, any number of which can be out of equilibrium. We note that advanced momentum-dependent methods have also been developed and used to treat neutrino oscillations in the early Universe [13][14][15][16].
We demonstrate our method in the context of the dark matter production near a sharp resonance, comparing our results with the ZOPLW approach and with the momentum-dependent method of [7], which uses a generalized relaxation time approximation for the numerically expensive backreaction terms. Our results validate the approximation scheme of [7] to its expected accuracy. Furthermore, we find that this scheme slightly overestimates the effect of elastic scattering channels, and this seems to be the case also with the truncated derivative methods used in [6,8]. We also apply our method in a sterile neutrino DM model including a singlet scalar mediator, first analyzed with simplified evolution equations in [10,11]. We find that neither including elastic interaction channels, nor adding new collision terms induced by a symmetry breaking changes the results appreciably. This verifies that the approximations used in [10,11] are robust and their results remain valid in the full solution.
The paper is organized as follows: In Sec. II we describe the discretization of the collision integrals. We then apply the developed methodology to two benchmark models. first, in Sec. III to the model where SM is extended with a real singlet scalar and then in Sec. IV to the model where the additional fields are a real singlet scalar and a sterile neutrino. In Sec. V we present our conclusions and outlook toward further work. Many details of the computations can be found in the Appendixes.

II. THE KINETIC EQUATION
In an expanding homogeneous and isotropic universe the Boltzmann equation can be written as where C coll are the collision terms describing the chemical and kinetic balances, H =ȧ/a is the Hubble parameter and a(t) is the scale factor. Expansion of the Universe is best quantified by integrating along the curves of constant comoving momentum, k = ap. Then, writing the Liouville operator in terms of k, the momentum derivative vanishes and we have where It is more natural to work with temperature instead of time, thus we define a dimensionless variable x ≡ m 0 /T , where m 0 is some reference scale, and T is photon temperature. To evaluate the Jacobian of this transformation, we use the adiabatic radiation era time-temperature relatioṅ s/s = −3H, where s = 2π 2 h eff T 3 /45 is the entropy density of the Universe and H = (4π 3 g eff /45) 1/2 T 2 /M P and g eff (T ) and h eff (T ) are the effective number of relativistic energy and entropy degrees of freedom. This implieṡ where We want to replace also the momentum with a dimensionless variable. Using again the adiabaticity condition, one finds (a/a 0 ) where k is the comoving momentum, p is the physical momentum, and a 0 is the scale factor evaluated at some reference temperature T 0 , which we set equal to the reference mass: T 0 ≡ m 0 . The Boltzmann equation in dimensionless variables becomes This Boltzmann equation, written in comoving variables, can be solved numerically by discretizing in variables x and ξ. Given such a discretization, the role of the parameters m 0 and T 0 is to tune the dimensionless variables to probe the desired temperatures and physical momenta. Before describing this process in detail, we must first carefully describe the structure of the collision terms C coll .

A. Collision integral
The collision term for generic two-particle interactions 12 ↔ 34 is given by where the integration measure over the phase space is We always assume that labels i = 1, 2, 3, 4 denote all internal degrees of freedom associated with a given distribution function with + (−) corresponding to the boson (fermion) case. Finally, |M 12→34 | 2 is the matrix element squared, summed or integrated over the internal degrees of freedom associated with the labels i = 2, 3, 4. The matrix element squared is also assumed to contain all relevant symmetry factors for the initial and final states. The collision integral C coll naturally splits into the backward and forward terms, given by and The phase space integration of the collision integrals has been studied in the context of neutrino astrophysics for massless neutrinos in [17] and later for nonzero neutrino masses in [18]. Similar methods were also developed, e.g., in [15,[19][20][21][22]. Here we follow the strategy of [18] to reduce the fully general, initially ninedimensional collision integrals down to four dimensions. The momentum dependence of the matrix elements prevents making further analytic simplifications.
Different from [18], we treat the forward and backward collision processes separately. This makes the numerical implementation more stable by avoiding the need to interpolate the unknown phase space distribution functions f (p, t) in between the integration grid points. Full details of the reduction are given in Appendix A. The final result for the reduced backward term (10) is given by (A11) and for the forward term (11) by (A15).

B. Discretization
We solve Eq. (6) numerically by discretizing the momentum grid ξ → ξ j , j = 1, . . . , N ξ , with uniform spacing in logarithmic scale. If the production processes spread over several orders of magnitude in temperature, this allows one to cover a sufficient range of momenta to reach the required accuracy. After discretization, the Boltzmann equation (6) becomes an initial value problem consisting of a coupled set of ordinary differential equations for f a (x, ξ j ) ≡ f aj (x), over some temperature range x, that must be solved simultaneously 1 for each degree of freedom a and the momentum mode ξ j , Here α(x) is the prefactor given in the rhs of Eq. (6), and the sum runs over all collision terms that contribute to evolution of f aj . Here we separated the degrees of freedom (labeled by a) from the discretized momentum variable. Indeed, each different particle species, and each helicity or polarization state within a species, in general has its own independent unknown distribution function, which the collision terms couple with each other. Some hierarchies between the interaction rates may allow simplifying the equation network, such as helicity equilibrium due to rapid helicity flips. This can be easily incorporated by imposing the degeneracies and introducing the corresponding averaged matrix elements. We will typically assume that initially f aj = 0 for the dark sector particle distributions. This is justified when we start early enough in time, i.e., high enough temperature, and it allows us to track to which degree each species thermalizes before it decays or its distribution freezes.
As an example, on collision term discretization we show how the backward term (10) is implemented. Although the Boltzmann equation (6) is solved in dimensionless momentum ξ, the collision term reduction in Appendix A is done in terms of the physical momentum p. The physical momentum p j corresponding to dimensionless momentum ξ j at a given temperature T = m 0 /x is then obtained by inverting Eq. (5), Then backward term (10) can be reduced to (A11), given in discretized form as where the discretized backward phase space factor is (15) and F 1jkl (x) is the angular integral over the matrix element squared defined in Eq. (A12). Here the superscripts 1, . . . , 4 denote the particle species (in the sense described above) involved in the 12 ↔ 34 process, f ai (x) is the value of distribution function of the particle "a" with momentum ξ i at temperature T = m 0 /x and energy E ai = (p 2 i + m 2 a ) 1/2 . The label j refers to the species 1, whose collision term we are computing and it is not summed over. The momentum space matrix structure f 2jkl of species 2 follows from the four-momentum conservation.
It is essential to note that the all matrices F 1jkl (x) can be precalculated and replaced by numerical fit functions for all relevant processes before solving the Boltzmann equations. This fitting procedure can be done very accurately and it is pivotal for the efficiency of the numerical code. The general flow of our implementation then is as follows: 1) Define theory, whose masses and couplings may depend on temperature.
2) Create grids and define the entropy and energy degrees of freedom functions h eff , g eff , g * . 6) Solve the Boltzmann system numerically using a suitable ordinary differential equation solver.
This formulation is generic enough to allow for dynamical changes that modify the parameters of the theory during the evolution, such as phase transitions. In the next sections, we show in detail the results of this implementation in simple hidden sector models connected with the SM via the Higgs portal.

III. FIRST BENCHMARK MODEL: THE SINGLET SCALAR EXTENSION
An extension of the SM by a real singlet scalar S, coupled with the SM Higgs doublet via the renormalizable operator |H| 2 S 2 , the "Higgs portal", provides a simple paradigm for a dark sector. Since its early introduction [23,24] this type of model building has started to gain more attention as benchmarks for experimental searches of particle dark matter [25][26][27][28]. Therefore, this model is a natural starting point for the tests of the computational method we have developed. Since the existing literature on this model is large and its phenomenology has been thoroughly exposed already, our discussion here will be brief; we will introduce only the necessary formulas and focus on the comparison of our approach with other approximate computation schemes. The singlet model is defined by the Lagrangian where the scalar potential is given by and the gauge interactions of the Higgs doublet H are contained in L SM . The stability of the potential requires that the quartic couplings λ S and λ H of the singlet and the Higgs fields are positive, but the portal coupling could be negative, as the stability of the potential requires just that λ HS > −2 √ λ S λ H . However, here we will only consider positive values of λ HS .
If the singlet scalar mass m s is just below half of the Higgs mass m H = 125.25 ± 0.17 GeV [29], the Higgs mediated inelastic processes are resonantly enhanced but elastic processes maintaining the kinetic equilibrium are not, so one would expect the nonequilibrium effects to be relevant. Indeed, if the singlet S is required to constitute all of the dark matter and remain compatible with the current direct detection experiments, its mass is constrained to m s ∈ [56, 62.5] GeV [28], where precision computations are required to address the dark matter phenomenology [7].
In our current, fully momentum-dependent setup, solving the singlet scalar relic density in this region entails solving the following Boltzmann equation: where s refers to the singlet, j = {τ, c, b, t, h, W, Z}, and f = {τ, c, b}. In this case we only need a dynamical equation for s, since all SM particles can be assumed to be in thermal equilibrium. To correctly include kinematics when the Higgs mediator in C I ss↔jj becomes on shell, we have included the on-shell Higgs decay C I h→ss as a separate contribution and take all Higgs mediators in C I ss↔jj to be off shell as described in Appendix C. Resonant inelastic processes deplete and overpopulate specific momentum states, which tends to bring the distribution function f s out of kinetic equilibrium. Elastic interactions, on the other hand, tend to restore the kinetic equilibrium, and if they are sufficiently fast, the standard thermal averaged treatment [4,5,28] suffices. However, the Higgs resonance is particularly sharp, 2 and one cannot a priori assume that the elastic processes can maintain the kinetic equilibrium to high accuracy. The issue has already been studied using momentumdependent methods [6,7], however, using some approximations in the treatment of the collision integrals.
In left panel of Fig. 1 we show the contours of Ω s h 2 = 0.1193 in the singlet mass and the portal coupling plane using different approximations. The results from the full computation implemented in this work are shown by the red crosses, while the standard thermal averaged result is shown by the solid blue line. The yellow circles correspond to the calculation using the momentumdependent, generalized relaxation approximation (GRA) method of Ref. [7], and finally, the calculation in Ref. [6], using truncated expansions for the elastic collision integrals, is shown by the black dots. GRA calculation is similar to the one described in this paper, except for the numerically expensive elastic backward term C sf ←sf , given by (10). In the GRA method, this is treated in a simplifying approximation; for more details, see [7].
When displayed in the logarithmic scale, all calculations appear to roughly agree. Plotting on linear scale (right panel in Fig. 1) reveals the significance of the deviations. The cause for the difference between the full and the GRA calculations is seen in Fig. 2, which shows the elastic collision integrals computed exactly (red solid line) and in the GRA (yellow dashed line). At high temperatures, the GRA method works well, but around the freeze-out temperature x ∼ 10 it starts to overestimate the elastic integral that enforces the kinetic equilibrium. Eventually the error becomes of order ∼ 2, but only well after the freeze-out x ∼ 40. This tendency was already noted in [7], and by construction the GRA scheme is not expected to work to a high precision for distributions that already are very close to thermal equilibrium. However, when one is close to equilibrium, the absolute magnitude of the error is already small, and GRA slightly improves on the thermal approximation. The results using a truncated expansion for the elastic collision integral from [6] are roughly comparable with the GRA.
Overall, we confirm that the effect of kinetic decoupling in the singlet model is not as dramatic as one might have initially guessed. All methods agree in the absence of resonant enhancement as they should. Even in the resonant region, by far the most important effect is to use the thermally averaged annihilation rate in the ZOPLW equations, first pointed out in [5]. The elastic interactions with the Standard Model particles are surprisingly efficient in keeping thermal equilibrium as observed in [7]. However, should a DM particle be identified in the resonance region, a precision calculation of its abundance  [7] (yellow circles), and complete momentumdependent computation with full backreaction terms (red crosses). Black dots are the full Boltzmann solution [30] corresponding to the method introduced in Ref. [6], for the case "QCD=A". Right: relative difference in Ωh 2 when λHS is fixed by the full calculation to yield Ωsh 2 = 0.1193 (red crosses on the left). Yellow circles correspond to relaxation approximation [7], blue squares to the momentum-independent method, and black dots correspond to the results of Ref. [6]. requires a full momentum-dependent calculation with an exact collision integral. Finally, let us note that the use of Maxwell-Boltzmann (MB) statistics in the inelastic collision integral causes about a 10% error [31] in final abundance. This is the case for all results in Fig. 1. We have checked that the corresponding error in elastic collision integrals only affects the final abundance by 1%.
To our knowledge, this is the first analysis of the singlet scalar in Higgs resonance with full elastic backreactions. This is undoubtedly due to the heavy numerical cost of the method. On average, we found that the computation times in this particular example (tested using a 6-core i7 laptop with 16 GB of RAM) scaled as follows in different approximations: thermal averaging runs took O(1) s, generalized relaxation approximation runs O(1−10) min, and the complete calculation with backreaction ∼ > 10 h per (m s , λ hs ) pair. Even faster and yet accurate methods of solving the ZOPLW equations exist [28] for the use of large scale parameter scans. Taking this hierarchy into account is obviously of paramount importance when making a choice of what method to use for a given problem.

IV. SECOND BENCHMARK MODEL: SINGLET SCALAR AND FERMION EXTENSION
As a second example, we consider an extension of the SM by a singlet scalar (S) and a singlet Dirac fermion (N ). The Lagrangian of the model is where we denote the SM Higgs doublet by H. Its gauge interactions are contained in the SM Lagrangian L SM , while the potential terms are contained in the extended scalar potential V (S, H) given by Eq. (17). In this model, the fermion N is a phenomenologically interesting candidate for cold DM [32,33]. We are interested especially in the keV mass range, where the nonequilibrium dynamics can be relevant [10,11], and the resulting nonthermal momentum distribution of DM may affect the formation of large scale structures. We focus on the question of whether highly nonthermal momentum distributions found in [10,11] survive when all elastic processes are included in the analysis. We will assume a mass hierarchy m N m S , with m N ∼ keV and m S ∼ O(10 − 1000) GeV. The N -fermion mass gets a contribution from nonzero vacuum expectation value (VEV) of the singlet scalar m N = µ N + y S . As the VEV can be quite large, we need to assume the Yukawa coupling to be tiny, y 1, to keep m N around keV scale. The vacuum structure is determined by the scalar sector of the theory. The field H is the usual weak doublet which has a VEV, denoted by v, along the neutral direction, φ 3 = v + φ. The VEV of the singlet field S is denoted by S ≡ w and we write S = w + σ. Inserting these parametrizations into Eq (19), setting the field fluctuations to zero and extremizing the full scalar potential leads to We use these conditions to eliminate µ 2 S and µ 2 H . This leads to the mass matrix for neutral scalars σ and φ, which is diagonalized by the transformation to the mass eigenbasis. We denote the mass eigenstates by h 1 and h 2 , so the explicit relation is We identify h 1 with the SM Higgs field and h 2 is a heavier scalar. Consistency with LHC data on Higgs couplings then requires sin θ ∼ < 0.23 [34,35]. We therefore set m 2 > m 1 = 125.25 GeV and consider the physical masses m 1 and m 2 to be input parameters. We then solve the couplings λ H and λ S and the mixing angle θ in terms of the physical masses,the vacuum expectation values and the portal coupling λ HS as Requiring sin(2θ) to be positive implies that 0 ≤ λ HS ≤ ≡ λ max HS . The Feynman rules following from the Lagrangian (19) are tabulated in Appendix B. In the special case of w = 0, the treatment is more straightforward, as the fields σ and φ are directly the mass eigenstates of the mass matrix. Without going into further details, we simply note that the Feynman rules of Appendix B can be directly applied also in this case by letting θ, w → 0. In the limit of vanishing Yukawa coupling and singlet scalar VEV, y, w → 0, the model reduces to the singlet scalar model from previous section.
With a slight abuse of notation, we denote the mass eigenstates by φ and σ, as this allows us to include the cases w = 0 and w = 0 simultaneously. Then we can summarize the above construction as follows: we have taken the masses m 1 ≡ m φ , m 2 ≡ m σ , the portal coupling λ HS , and the vacuum expectation values w and v as the input parameters, and express other Lagrangian parameters in the scalar sector in terms of these. Furthermore, we fix v = 246 GeV. Thus, the free parameters in this theory are {m σ , λ HS , w, y, m N }.

A. DM production processes and coupled Boltzmann system
The Lagrangian (19) allows for various production processes for the N and σ. Processes of order O(y 2 ) are negligible and the relevant contributions under our assumptions are summarized in Table I. Because of nonzero Yukawa coupling and assumed mass hierarchy, eventually all produced σ scalars will decay into N fermions, which remains as a stable relic. Production of σ scalars is determined by the portal coupling λ HS . Direct production of N fermions from a SM heat bath is allowed by a nonzero mixing angle sin θ between the scalars, but remains subdominant for the allowed small mixing angles. Therefore, the production of N fermions proceeds mostly via σ-scalar decays, which is itself produced from a SM heat bath and whose number density can freeze (either via freeze-in or freeze-out mechanism) before it fully decays.
To obtain the momentum distribution function for σ scalar and N fermion we must solve the following set of

Always open
Open  coupled Boltzmann equations:  [36,37]. Here we are following the treatment of [38,39]; see Appendix C for more details and discussion.
Different from previous treatments, we have also accounted for the three-and four-body final states from virtual boson decays using methods described in [28], as well as the one-loop corrections for quarks in the σσ ↔ jj channel. Accounting for virtual boson decays and QCD one-loop corrections describe the SM states more accurately and slightly increase the C I σσ↔jj contributions in Eq. (29). This is good to keep in mind when comparing our results to, e.g., Ref. [11], as in the case of σ freezing out this slight increase causes the σ to follow the SM heat bath a bit longer and slightly suppresses the final fermion distribution.

B. Results and discussion when w = 0
We first set the VEV of the singlet scalar to zero, so that the scalars do not mix. This leaves us with processes on the left column of Table I. This setting is equivalent to the one studied in Ref. [11], except that we have included  [11] (black crosses). Note that our convention to λHS differs from [11] by a factor of 4. Darker colors refer to later times and the black solid curve is the final frozen-in form. The low-end momentum tail is produced at temperatures much higher than the electroweak phase transition temperature, T TEWPT = 150 GeV. The temperature grid range is log 10 (x) ∈ [−3.5, 3] and we used m0 = T0 = m h as the reference scales in (5) and (6).
the elastic processes σf ↔ σf , σσ ↔ σσ, and N σ ↔ N φ, which tend to suppress the nonthermal component in the momentum distribution of N fermion. The two first processes can also lengthen the freeze-out time of σ field, thus allowing it to be Boltzmann suppressed more before it freezes out and decays, which can reduce the late time production of N fermions. The two-peaked nonthermal momentum distribution found in [11] results from N being produced at two separate temperature scales (see [12] for a comprehensive study). Hence, reducing the production at either temperature scale could prevent the momentum distribution from forming the double peak structure. The last two processes tend to restore the kinetic equilibrium by reducing the nonthermal component momentum distribution. In practice, we find their effect to be negligible. Figure 3 shows our solution for the momentum distribution function of the fermion N for parameter values {m σ = 60 GeV, m N = 20 keV, λ HS = 4 × 10 −6.86 , y = 10 −8.69 , w = 0}. These values are chosen to correspond to the case presented on the right-hand panels of Fig. 5 in Ref. [11] and we have included their final result as black crosses in our Fig. 3. Clearly the elastic processes are not sufficiently fast to suppress the formation of the nonthermal momentum distribution. The reason is quite clear a posteriori: The second peak at high momenta forms only when σ production from the SM heat bath has already stopped and the remaining scalars decay into N fermions. At this point, the SM Higgs is also heavily Boltzmann suppressed. The elastic channel N σ ↔ N φ, which would most effectively thermalize the momentum distribution of N fermions, is therefore suppressed due to absence of both scalars in the thermal bath.
We have also checked the other light (m σ < m h /2) and heavy (m σ > m h /2) scalar (corresponding to feebly and weakly interacting massive particles) scenarios as discussed in [11] and found that our results agree with theirs to the same extent as in the above example. In summary, we find that the model presented in [11] is indeed inherently nonthermal in parts of its parameter space when assumptions y 1 and w = 0 hold. Our numerical results for the distribution f N (ξ, x) differ noticeably from [11], but the difference does not come from elastic interactions, but from our more accurate evaluation of the inelastic interaction rates.
C. Results and discussion when w > 0 We now let the singlet obtain a nonzero VEV w > 0. The nonzero mixing angle induced by the broken symmetry is given by Eq. (27). It increases the σ and N particle production rates and opens new elastic channels for the N field. All relevant reaction channels are given in Table I, where also the channels on the right column now contribute.
In Fig. 4 we compare our results in the w > 0 case with the previous w = 0 case (shown with the green curve) to see the effect of additional channels and the increased interaction rates on the nonthermal double peak structure of the momentum distribution of N fermions. We have chosen conservative values for the singlet VEV: w = 2 × 10 3 GeV (blue curve) and w = 2 × 10 6 GeV (red curve). In both, cases we find a strong suppression on the amplitude of the nonthermal second peak, which strongly increases for increasing w. However, this is again not due to elastic processes being effective in redistributing the quanta, but due to the fact that the inelastic interactions between the σ scalar and SM heat bath strongly increase for increasing values of w, as illustrated in Fig. 5. As a result, the number density of σ particles becomes more and more Boltzmann suppressed before they finally decay to N fermions. This then suppresses or removes the high-momentum out-ofequilibrium tail from the N distribution. The mechanism is clearly illustrated in Fig. 6, where we show the integrated yields of the singlet scalar fields (dotted lines) as a function of m 0 /T for the same choices of parameters as in Fig. 4. The effect on integrated neutrino distributions (dashed lines) is much smaller than to their momentum dependence, but it shows the right tendency as function of the σ-field abundance.
We found that the elastic interactions had negligible effect on the N distribution in comparison with the Boltzmann suppression discussed above. Increasing the VEV of the singlet does result in stronger elastic rates for the fermion N , but they still fail to restore the thermal equilibrium. This is partly because the having w > 0 also increases the overall production of N fermions via inelastic and decay channels, which is the root cause for the nonthermal distributions. We then find that even in the  Fig. 3, and the blue dotted line shows an example of forced thermalization when the elastic rates are increased by hand, which we show only as a proof of concept to show that the kinematics work as indented and the system then tends toward thermal distribution.
w > 0 case the elastic processes are never sufficiently fast to thermalize the fermion momentum distribution and the model (19) is inherently nonthermal when y 1 and m N m S . However, if a phase transition takes place in the singlet sector before the singlet field freezes out, the predictions for the amplitude of the out-of-equilibrium component in N distribution change dramatically for a given set of parameters, in comparison to the w = 0 case studied in [11]. However, this change is not due to increased elastic interactions, but due to change in inelastic rates, and to discover this effect the high-accuracy Boltzmann codes developed in this work are, in fact not necessary.
The inefficacy of the elastic rates in the examples studied above results from the smallness of the adopted values for the y coupling and this does not imply that elastic interactions were irrelevant in general. Instead of constructing a model just to make the elastic rates important, as a proof of concept, we simply increase the elastic process N f ↔ N f by hand to see how large it must be for a significant thermalization by elastic channels to take place. For {m σ = 60 GeV, m N = 20 keV, λ HS = 4 × 10 −6.86 , y = 10 −8.69 , w = 2 × 10 3 } we scaled the matrix element M N f →N f → 10 11 ×M N f →N f . The nonscaled result is shown as a solid blue curve and the result after scaling as a dotted blue curve in Fig. 4. Thus, in this case, at the level of cross sections, the elastic rates are roughly ∼ 10 −22 times too small to significantly thermalize the system.
Let us finally note that, while the introduction of the phase transition changes the predictions significantly for a given portal coupling, the same out-of-equilibrium distribution can be obtained in the spontaneously broken case for a different portal coupling. That is, there are degenerate subspaces in the (λ HS , y; w) parameter space, where the same nonthermal momentum distribution can be obtained. In particular for a fixed Yukawa y, essentially only the interaction rate between the singlet scalar and SM heat bath is important in determining the degree to which the momentum distribution becomes nonthermal. As this rate is determined by a combination of λ HS and w, we can always find such a (w , λ HS ) pair that the original nonthermal behavior for (λ HS , y; w = 0) parameters is obtained with (λ HS , y; w > 0). This is shown in Fig. 7, where we plot the distributions with y = 10 −8.69 for the cases {w = 2 × 10 3 GeV, λ HS = 4 × 10 −6.86 } and {w = 2 × 10 6 GeV, λ HS = 4 × 10 −8 }, respectively. In each case we find exactly the same momentum distribution, apart from small deviations in the very small momenta.

V. CONCLUSIONS
We have presented a computational method that is generally applicable for solving the coupled set of Boltzmann equations for phase space distribution functions of cosmic relics. Similar techniques have been developed earlier in the context of the neutrino astrophysics [13][14][15][16], but until now they have not been carefully adapted to be used in the dark matter abundance calculations. For earlier implementations that use different levels of approximations for the collision integrals, see [6,7]. One of the main results is the compilation of these methods into a simple and generic form that provides a suitable starting point for their wider utilization in cases, where the standard thermal averaging methods fail and more accurate predictions in DM models are needed.
We demonstrated our method in two models known to be sensitive to the momentum dependency of the phase space distribution. Our first example concerned the freeze-out of a singlet scalar dark matter coupled to SM via the Higgs portal near Higgs resonance. In this case, the kinetic equilibrium approximation required by the usual ZOPLW equation cannot be assumed a priori. We performed a full momentum calculations, comparing our results to the ZOPLW approximation and two earlier momentum-dependent calculations, where further approximations of the form of the elastic collision integrals were made.
We found that the approximation methods of Refs. [6,7] are in good agreement with our full results. Both methods present some improvement over the ZOPLW results even very near equilibrium, although they both slightly (the former a little more) overestimate the elastic rates.
As another example, we considered the model presented in [10,11]. Using our methodology we extended their calculations to include the elastic processes and possible mixing between the two scalar states as a result of a phase transition in the singlet sector. We found that this could significantly alter the predicted size of the nonequilibrium distribution, the more so the larger the VEV of the singlet field. However, this difference was mainly due to changes in the inelastic rates, caused by the phase transition. Elastic rates turned out to be inefficient and even when they were included, the momentum distribution of the singlet fermion remains inherently nonthermal. Overall, we find that, in a vast majority of cases, the momentum averaged methods work surprisingly well. ACKNOWLEDGMENTS We thank T. Bringmann for correspondence and sharing data related to Fig. 1

NOTE ADDED
Recently, a paper presenting similar methods, also based on [18], appeared in [40]. Our results agree qualitatively with theirs.
, and the function F (p 1 , p 3 , p 4 ) contains the squared matrix element integrated over the angles and fixing kinematics, where we set z = cos α and z ± ≡ (−b ± √ b 2 − 4ac)/2a. This function contains all process specific dynamical information and since it is independent of the distribution functions, it can be computed and fitted before solving the Boltzmann equations. In general, the matrix element squared is a function of s = (p 3 + p 4 ) 2 and t = (p 1 − p 3 ) 2 , which depend on the angles and momenta as follows: In the particular case where the matrix element in (A12) does not depend on cos α (a pure t-channel process), the d cos integral can be reduced to a one-dimensional integral, (A14) Forward term. In the forward term, given by Eq. (11), we integrate first over d 3 p 4 , which leaves us with the delta function δ(p 2 4 − m 2 4 ). Paralleling the backward term reduction, we eventually obtain where f 4 is evaluated at E 4 = E 1 + E 2 − E 3 and F (p 1 , p 2 , p 3 ) has an identical expression to the righthand side of (A12), where one replaces everywhere where Q = m 2 1 + m 2 2 + m 2 3 − m 2 4 and, moreover, γ = E 1 E 2 −E 1 E 3 −E 2 E 3 and = p 1 p 3 cos θ and κ = p 2 1 +p 2 3 . Now the Mandelstam variables must be written as s = (p 1 + p 2 ) 2 and t = (p 1 − p 3 ) 2 so that If the matrix element is again independent of cos α, the result (A14) applies also as such, after replacing (a, b, c) → (a , b , c ). Note that the forward and backward collision integrals (A11) and (A15) are valid for general Bose-Einstein and Fermi-Dirac statistics.

Special case: Inelastic 2-2 scattering toward/from equilibrium in MB statistics
As stated in the main text, we are using the MB statistics throughout. This is not needed for our computation of the elastic rates, but to keep the computation time associated with the large number of inelastic interactions with the SM states manageable. In this case, the final states are in equilibrium, and we can reduce the ninedimensional integral down to one-dimension. Overall, using the MB statistics amounts to about 10% error in the overall magnitude of the elastic collision integral [31], which should have but a very small effect on the final abundance. Indeed, we checked that scaling the elastic collision terms by a factor 0.9-1.1 caused only a 0.7% change in the final abundance. Then, working under the assumption f 1 and enforcing the detailed balance, we can write the phase space factor (9) as The collision can then be written as where n runs over different equilibrium states, σ where λ(a, b, c) = (a − b − c) 2 − 4bc is the Källén kinetic function. As explained in [7], one can reduce the integral over d 3p 2 to a single integral over s, with s ± = m 2 1 + m 2 2 + 2E 1 E 2 ± 2p 1 p 2 . Again, this function can be evaluated and fitted for each collision channel before one attempts to solve the dynamical Boltzmann equations, which gives a dramatic boost in numerical efficiency.

1-2 decays and fusions
A similar reduction that was carried out above for the 2 ↔ 2 scatterings, can be performed for the 1 ↔ 2 processes. We shall assume that either the decaying particle or the decay products are in thermal equilibrium.
Decay from (fusion to) equilibrium f eq A ↔ f 1 f 2 . Assume we are tracking the species labeled as "1", while the "2" species is arbitrary and "A" species follows the equilibrium. Then, where the phase space integration is denoted by and the distribution factor Λ is given by where we again assumed that f 1 and applied the detailed balance for the equilibrium state. A similar procedure as in the previous section eventually gives C eq A−12 (p 1 , t) = where we used the fact that the matrix element for the decay process is a constant and defined with In this case one only needs to compute the matrix element as a function of the masses of particles involved.
Decay to (fusion from) equilibrium f 1 ↔ f eq A f eq B . Now assume we track the species "1" while the arbitrary species "A" follows the equilibrium. After similar steps as above, we get an even simpler expression C eq 1→AB (p 1 , t) = where v = λ 1/2 (m 2 1 , m 2 A , m 2 B )/m 2 1 and where we used the fact that the matrix element squared (in tree level) is always a constant.  [39]. The split is done for the propagator and the resulting difference shown for the propagator squared. Full propagator D(p 2 ) is shown with the black solid line, the off-shell part DH(p 2 ) by the dashdotted line, and the on-shell part A(p 2 ) by dash-dotted lines. The vertical lines mark the resonance width p = √ m 2 ± mΓ, which coincides with the maxima of the off-shell propagator.
One can argue for this prescription by noting that, for a simple matrix element with no mixing between different channels, |M| 2 ∼ D(p 2 )D * (p 2 ) = D 2 H (p 2 ) + A 2 (p 2 ). (C4) To the lowest order in a small but finite Γ, the square of the spectral function can be replaced by The prescription (C3) is thus approximatively the same as (C2). However, this argument fails in the presence of interference terms and the prescription (C3) has been found to give superficial negative cross sections [38], which never happens in the subtraction scheme (C2). In Fig. 8 we show the split of the square of the propagator function into the on-and off-shell contributions according to (C2) in a representative case.