{'fraction_non_alphanumeric': 0.03541315345699832, 'fraction_numerical': 0.009724564362001125, 'mean_word_length': 4.376397703233605, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Thermal Abundance of Semi-Relativistic Relics

Approximate analytical solutions of the Boltzmann equation for particles that are either extremely relativistic or non-relativistic when they decouple from the thermal bath are well established. However, no analytical formula for the relic density of particles that are semi-relativistic at decoupling is yet known. We propose a new ansatz for the thermal average of the annihilation cross sections for such particles, and find a semi-analytical treatment for calculating their relic densities. As examples, we consider Majorana- and Dirac-type neutrinos. We show that such semi-relativistic relics cannot be good cold Dark Matter candidates. However, late decays of meta-stable semi-relativistic relics might have released a large amount of entropy, thereby diluting the density of other, unwanted relics.


Introduction
The accurate determination of cosmological parameters by up-to-date observations, most notably by the Wilkinson Microwave Anisotropy Probe (WMAP) [1], increases the importance of quantitative predictions. In particular, the estimate of the cosmological relic abundances of particle species is essential, for the history of the universe depends on these quantities. One of the most important examples is the cosmological abundance of dark matter [2,3], whose mass density is found to be [1] Ω DM h 2 = 0.1099 ± 0.0062 , where h ≃ 0.7 is the scaled Hubble parameter in units of 100 km sec −1 Mpc −1 and Ω is the mass density in units of the critical density.
The abundance of some particle species is determined by solving the Boltzmann equation, which describes the change of the particle number caused by particle reactions as well as by the expansion of the universe [2]. However, there is no analytical general solution of this nonlinear differential equation, and therefore one needs to solve the equation numerically in many cases. In early studies, approximate analytical formulae have been found for the relativistic [4,2] and non-relativistic [5,6,7] regimes.
Stable or long-lived weakly interacting massive particles (WIMPs) with weakscale masses are examples of cold relic particles, which decouple from thermal equilibrium when they are non-relativistic. In standard cosmology, decoupling of WIMPs occurs in the radiation-dominated (RD) era after inflation, and analytic approximate formulae for the WIMP relic abundance have been derived [5,6,7]. In the opposite limit, where decoupling of particles from the thermal background occurs when they are relativistic, the relic abundance is approximated by its equilibrium value at the decoupling temperature, and not sensitive to details of its freeze-out [2]. The resulting relic density can therefore easily be computed analytically. On the other hand, no analytical treatment to calculate the relic abundance in the intermediate regime is known yet.
In this paper, we revisit the relic density of particles χ that decouple from the thermal bath when they were semi-relativistic, i.e. at freeze-out temperature T F ∼ m χ . Assuming that the Maxwell-Boltzmann distribution can be used for all participating particles, we introduce an expression for the thermal average of the χ annihilation cross section which smoothly interpolates between the extremely relativistic and non-relativistic regimes. It is shown that our new ansatz is capable of reproducing the exact thermally-averaged annihilation cross section with accuracy of a few percent. Given this approximated cross section, we can define the freeze-out temperature by comparing the annihilation rate to the expansion rate. The assumption that the comoving χ density remains constant after freeze-out turns out to be a good approximation for the relic abundance of semi-relativistic particles.
We also discuss the roles such semi-relativistic particles could play in realistic cosmological scenarios. It should be emphasized that the abundance of semirelativistically decoupled relics tends to be large because it is only very mildly Boltzmann suppressed. We point out that scenarios where semi-relativistically decoupling particles form the Dark Matter have problems with structure formation, Big Bang Nucleosynthesis (BBN) and/or laboratory measurements. Nevertheless, such semirelativistic relics can be useful for diluting the density of other, unwanted relics by late-time out-of-equilibrium decay [8,9]. As an example, we investigate a scenario of decaying sterile neutrino that is assumed to depart from thermal equilibrium when it is semi-relativistic, in sharp contrast to non-thermal sterile neutrino scenarios [10]. Thermal equilibrium is attained by introducing some higher-dimensional operator. It is illustrated that an enormous amount of entropy can be produced without spoiling the successful BBN prediction of the light element abundances.
This paper is organized as follows: In Sec. 2 we begin by reviewing briefly the method to calculate relic densities for relativistic and non-relativistic particles. In Sec. 3 we explain the new formalism which is applicable for all freeze-out temperatures in case of S-and P -wave cross sections. The way to calculate the freeze-out temperature is also shown. In Sec. 4, the possibility for semi-relativistic particles to have the observed dark matter relic density of Ω DM h 2 ≃ 0.1 is considered. Then, we discuss the amount of entropy produced by the decay of unstable semi-relativistic species that decay in less than a second. Finally, Sec. 6 is devoted to summary and conclusions. Some properties of modified Bessel functions are described in Appendix A. In Appendix B we argue that the use of the Maxwell-Boltzmann distribution in the definition of the thermally averaged cross section only leads to a small mistake in the final relic density.

Relic abundances in the non-relativistic and relativistic limits
In this Section we briefly review the standard analytical approximations for evaluating the relic abundance of hypothetical particles χ [2,5,7]. These are applicable to particles that were either non-relativistic or extremely relativistic at freeze-out.
The number density n χ is obtained by solving the corresponding Boltzmann equation [2]. For the moment, we assume that single χ production and χ decay are forbidden by some symmetry or adequately suppressed.
where n χ,eq is the χ equilibrium number density, σv is the thermal average of the annihilation cross section multiplied by the relative velocity between the two annihilating χ particles, and H is the Hubble expansion rate of the universe. The second term on the left-hand side describes the dilution caused by the expansion of the universe; the first (second) term on the right-hand side decreases (increases) the number density due to annihilation into (production from) other particles, which are assumed to be in complete thermal equilibrium.
It is useful to express the above Boltzmann equation in terms of the dimensionless quantities Y χ = n χ /s and Y χ,eq = n χ,eq /s. The entropy density is given by s = (2π 2 /45)g * T 3 , with g * being the effective number of relativistic degrees of freedom and T the temperature of the universe. We also introduce the dimensionless ratio of χ mass m χ to the temperature, x = m χ /T . Assuming an adiabatic expansion of the universe, Eq.(2) can then be rewritten as [2] dY The generic picture of χ decoupling from the thermal bath is as follows. After inflation, the universe becomes radiation-dominated with expansion rate where M Pl = 2.4 × 10 18 GeV is the reduced Planck mass. The reheat temperate is assumed to be high enough for χ particles to reach full thermal (chemical as well as kinetic) equilibrium. * Thermal equilibrium is maintained as long as the interaction rate Γ = n χ σv is larger than the Hubble expansion rate H. As the temperature decreases, the interaction rate decreases more rapidly than the expansion rate does.
When the interaction rate falls below the expansion rate, χ is no longer kept in thermal equilibrium and the comoving number density Y χ becomes essentially constant.
This transition temperature is referred to as the freeze-out temperature T F .
Analytical expressions for the resulting χ relic density are known for the cases where decoupling happens when χ is non-relativistic ( . We discuss these two limiting cases in the following Subsections.

Relativistic case
First, consider the case where particles χ decouple when they are ultra-relativistic (x F ≪ 3). In this case the equilibrium number density to entropy ratio Y χ,eq (x) depends on the temperature only through the number of degrees of freedom g * of the thermal bath. Therefore, the final relic abundance is to very good approximation equal to its equilibrium value at the time of decoupling: where g eff = g χ (for bosons) , with g χ being the number of internal (e.g., spin or color) degrees of freedom of χ. Following the conventional notation, we express the χ relic density as Ω χ = Pl is the present critical density of the universe, and s 0 ≃ 2900 cm −3 is the present entropy density. This yields It should be noted that the relic density is simply proportional to the mass of the particle in the relativistic case. * The case where the reheat temperature is too low for χ to attain chemical equilibrium has been discussed in [11,12].

Non-relativistic case
For the non-relativistic case where x F ≫ 3, the relic abundance strongly depends on the freeze-out temperature T F because the equilibrium abundance Y χ,eq (x) is exponentially suppressed as the temperature decreases. The temperature dependence of the thermal average of the annihilation cross section is obtained using the Taylor expansion in powers of the velocity squared: The numerically-evaluated correct relic abundance is reproduced with an accuracy of a few percent using the approximate analytic formula For WIMPs with electroweak scale mass, freeze-out occurs at x F ≃ 20. The corresponding scaled relic density is then given by Note that the relic density of a non-relativistic particle is inversely proportional to its annihilation cross section, but does not depend explicitly on its mass.

Abundance of semi-relativistically decoupling particles
In the previous Section, we reviewed the known relativistic and non-relativistic approximate formulae for the relic abundance. The main aim of this Section is to find a simple method applicable between the two regimes: an analytic estimate of the relic abundance of semi-relativistically decoupling particles (x F ∼ 3).
One of the key quantities that determine the freeze-out temperature is the thermally-averaged cross section. In the non-relativistic case, expressions for the equilibrium number density as well as for the thermally-averaged cross section times velocity are rather simple. For a semi-relativistic particle, however, the thermallyaveraged cross section involves multiple integrals and cannot be expanded with respect to the velocity nor to the mass. Here we discuss a method of approximating the thermally-averaged cross section, by interpolating between its relativistically and non-relativistically expanded expressions. We employ the Maxwell-Boltzmann distribution for the equilibrium number density [7]: The thermal average of the cross section is then obtained as [7] σv = 1 8m 4 where K 1 (x) and K 2 (x) is the first and second modified Bessel function of the second kind; some properties of these functions are given in Appendix A.
At first glance the use of the Maxwell-Boltzmann distribution seems improper for particles that are semi-relativistic at decoupling, let alone for ultra-relativistic particles. However, we will argue in Appendix B that this should still yield accurate results for the final relic density. This is partly due to cancellations between the numerator and denominator of Eq.(12), and partly due to the fact that the final result for Ω χ becomes less sensitive to x F as x F decreases.
As examples of the annihilation of particles, we consider the pair annihilation processes of Dirac and Majorana fermions (e.g. neutrinos) into a pair of massless fermions, χχ → ff [4,13,14]. It should, however, be noticed that this assumption includes more general cases of any other species annihilating from S− or P −wave initial states. In a renormalizable model, the annihilation is mediated by some heavy particle: for example, a Z-boson with tiny coupling with χ, or a new spin-1 boson U [15]. We assume that the mass of this exchange particle is much larger than m χ , so that the annihilation amplitude can be described through an effective four-fermion interaction.
The annihilation of two Dirac fermions proceeds from an S−wave, and the resulting cross section can be parameterized as where s is the center-of-mass energy squared. G denotes the coupling constant of the four-fermion interaction (e.g. the Fermi coupling constant, G F = 1.17 × 10 −5 where p A,B and E A,B are the four-momenta and energies of the two incident particles labeled A and B. The resulting thermally averaged cross section is given by Its relativistic and non-relativistic limits read respectively. A general expression for σ S v should reproduce these results for x → 0 and x → ∞, respectively. A simple possibility is It should be noticed that this choice is not unique. Let us turn to the case of the annihilation from a P −wave initial state, which is e.g. true if χ is a Majorana fermion. Eq.(13) should then be replaced by Thermal averaging leads to Following the same steps as in the S−wave case, we find the following interpolation: Figure 1 shows the ratio of the approximate to the exact cross section σv app / σv for the S− (solid line) and P −wave (dashed) cases. Note that this ratio depends only on x. We see that our approximate expressions reproduce the exact ones with accuracy of better than 2% (0.5%) for annihilation from an S− (P −)wave, even in the semi-relativistic region (x ∼ 3).
Using Eqs. (17) or (20) instead of the exact expressions (15) or (19)   We define the freeze-out temperature by equalizing the interaction rate and the We then simply assume that Y χ does not change after decoupling from the thermal † This definition of the freeze-out temperature should not be used for the prediction of the relic abundance through Eqs. (9) and (10).  bath, so that Let us see to what extent this method can reproduce the correct relic abundance.
As an example, we consider the pair annihilation of neutrino-like particles via the mediation of the weak SM gauge bosons ‡ . In the left frame of Fig. 2 we plot the relic abundance Ω χ h 2 as function of m χ . The solid curves show predictions for the relic abundance obtained by solving the Boltzmann equation (3) numerically, while the dashed curves have been obtained using the analytic approximation described above. The upper curves are for Majorana fermions annihilating from a P −wave, and the lower curves are for Dirac fermions case annihilating from an S−wave. We take G = G F = 1.17 × 10 −5 GeV 2 , g * = 10 and g χ = 2 § . The right frame shows the corresponding values of x F . ‡ Our primary concern here is to test our approximation for the relic abundance, and thus, as an illustration, we take such an unrealistic setup. § Since we concentrate on the ratio of the exact to the approximate relic densities, we discard the temperature dependence of g * , which is basically an overall factor.
This figure shows that our very simple analytical treatment reproduces the correct relic density with an error of at most 20% (5%) for semi-relativistically decoupling particles annihilating from an S− (P −)wave. Not surprisingly, our treatment becomes exact for particles that are relativistic at decoupling. ¶ In the Dirac case, our approximation coincides with the exact relic abundance for m χ = 1 GeV, corresponding to x F ≃ 16; however, the deviation becomes larger again for larger m χ .
Therefore, our method is not applicable for the entire region of cold relics. Instead, one could switch to the usual non-relativistic treatment described in Sec. 4 Semi-relativistic dark matter?
As a first application, let us analyze whether semi-relativistically decoupled particles (x F ≃ 3) can be a dark matter candidate, whose cosmological abundance should be Ω χ h 2 < ∼ 0.1. The final number density of such a particle is of order of Y χ,eq (x ≃ 3) ∼ 10 −2 . Combining this value with the observed amount of dark matter, the upper bound of the semi-relativistic particle turns out to be m χ < ∼ 100 eV, which would thus decouple at temperatures of a few dozen eV. These particles would therefore still be ultra-relativistic during the formation of 4 He. Moreover, the effective coupling would have to be very large, G ∼ 10 3 GeV −2 . Such a scenario is therefore tightly constrained.
Note that a χ particle with m χ ∼ 100 eV could only annihilate into light neutrinos. * Moreover, the exchange particle would also need to be quite light, with mass < ∼ 30 MeV if all couplings are ≤ 1.
The simplest case is that of a real scalar χ. Since it only adds a single degree of freedom, its presence would not be in serious conflict with current BBN constraints ¶ We will see in Appendix B that one should use the correct Fermi-Dirac or Bose-Einstein expression for n χ,eq in Eq. (22) in order to accurately predict the relic density for x F < ∼ 1. * The neutrinos themselves aren't in thermal equilibrium with the photons any more at T < ∼ 100 eV. However, as long as T ≫ m ν they would still have a thermal distribution. The Boltzmann equation (2) should therefore still be applicable, if T is taken to be the temperature of the neutrinos, which is somewhat lower than the photon temperature. [16]. In a renormalizable theory, χχ → νν could then proceed either through exchange of a fermion in the t− or u−channel, or through boson exchange in the s−channel.
In the former case, the exchange particle would have to be an SU(2) doublet, if the low-energy theory only contains left-handed, SU(2) doublet neutrinos. The presence of such a light SU(2) doublet fermion is excluded by LEP data. In principle the light neutrinos might also be Dirac particles, allowing χχ → ν R ν R annihilation via exchange of a singlet fermion, possibly ν R itself. However, one would then need ν R to also be in thermal equilibrium, increasing the number of additional degrees of freedom present during BBN to an unacceptable level.
For a real scalar χ, s−channel exchange could only proceed through another scalar φ. However, a scalar φ can only couple to ν L ν R or its hermitean conjugate.
This scenario would therefore again require ν R to have been in equilibrium. We therefore conclude that such a light χ particle cannot be a real scalar. † If χ is a complex field, one needs g χ = 2, which is only marginally compatible with BBN [16]. In principle, χχ → νν could then proceed through t− or u−channel exchange of an SU(2) singlet. However, then either χ or this light exchange particle would have to carry hypercharge, so that it would have been produced copiously in Z decays. The argument against s−channel exchange of a scalar is the same as for real χ.
However, a complex χ, either a scalar or a Weyl fermion ‡ , could couple to a new gauge boson U, which in turn could couple to ν L ν L . This new boson would contribute three additional bosonic degrees of freedom at T > ∼ m U . Since we already added two degrees of freedom in χ, consistency with BBN would require m U > ∼ 1 MeV. This in turn would require the coupling of U to left-handed neutrinos to exceed 0.01, assuming its coupling to χ is perturbative. By SU(2) invariance, U would have to couple with equal strength to left-handed charged lepton. This combination of U boson mass and coupling is excluded, by a large margin, for both the electron and muon family by measurements of the respective magnetic moments [15]. No analogous measurement exists for the third generation, so a U−boson with few MeV mass coupling exclusively to third generation leptons and χ particles might still be † Of course, this argument does not exclude the possibility that a much heavier real singlet scalar χ could be cold Dark Matter [17].
‡ A massive two-component Weyl fermion can equivalently be described by a four-component Majorana fermion.
compatible with laboratory data. However, given that µ and τ neutrinos are known to mix strongly [18], a gauge invariant model where U couples to ν τ but does not couple to muons is difficult, if not impossible, to construct.
Finally, such a light χ particle would form hot, or at least warm, Dark Matter.
This possibility is strongly constrained by observations of early structures in the universe, in particular the "Lyman−α forest" [19]. Such a χ particle could thus at best form a sub-dominant component of the total Dark Matter.
In combination, these arguments strongly indicate that semi-relativistically decoupling particles should not be absolutely stable. In the next Section we will show that such particles may nevertheless have a role to play in the history of the Universe, if they are metastable.

Entropy production by decaying particles
In this section we demonstrate that semi-relativistically decoupling particles can be useful for producing a large amount of entropy, which could dilute the density of other relics to an acceptable level. Examples of such relics are decaying gravitinos, which can lead to problems with Big Bang nucleosynthesis [20], or supersymmetric neutralinos, whose relic density often exceeds the required Dark Matter density by one or two orders of magnitude [21]. The density of such relics will be diluted only if the entropy is released after they decouple from the thermal bath. This will simultaneously dilute any pre-existing baryon asymmetry. One thus either has to increase the efficiency of early baryogenesis, or introduce late baryogenesis after the release of the additional entropy. Both possibilities can be realized in the framework of Affleck-Dine baryogenesis [22].
Generally [8,9], out-of-equilibrium decays of long-lived particles can only produce a significant amount of entropy if the decaying particle dominates the energy density of the Universe prior to its decay. The abundance of non-relativistically decoupling particles is suppressed by a factor e −x F , hence their contribution to the energy density is small at decoupling. However, after decoupling their energy density only drops like R −3 ∝ T 3 , while that of the dominant radiation component decreases like T 4 as the Universe cools off. Thermally produced particles can therefore dominate the energy density of the Universe only at temperature T ≪ e −x F T F .
Significant entropy production by the late decay of nonrelativistically decoupling particles is therefore only possible if they are simultaneously very massive and quite long-lived. For semi-relativistic particles, on the other hand, the abundance at decoupling is large and thus a significant amount of entropy can be produced even if their mass is small, since their density will become dominant quite soon after decoupling.
Let us consider the out-of-equilibrium decay of long-lived particles which semirelativistically decoupled from the thermal background. For simplicity we work in the instantaneous decay approximation, i.e. we assume that all χ particles decay at time t d = τ χ , where τ χ is the lifetime of χ. While this approximation does not describe the time dependence of the entropy (or temperature) for t ∼ τ χ very well, it does reproduce the entropy enhancement factor, i.e. the entropy at t ≫ τ χ , quite accurately. We assume that χ particles were in full thermal equilibrium for sufficiently high temperatures in the RD epoch. When the temperature decreased to T = T F ≃ m χ , the χ number density n χ froze out. At decoupling, χ particles contributed a few percent to the total energy density of the universe; however, as noted earlier, the ratio of the radiation and χ energy densities decreased by a factor T d /T F = t F /τ χ between decoupling and decay of χ; here T d refers to the temperature at time t = τ χ , just prior to χ decay. If τ χ ≫ t F , the χ energy density at the time of the χ decay is well approximated by ρ χ,d = m χ n χ,d , and dominated over the radiation. In this case, the ratio of the final entropy density s f after the decay to the initial entropy density s i before the decay is given by [8] for s i ≪ s f . Here Y χ,d = n χ,d /s i is proportional to the χ abundance just prior to its decay.
In the light of the BBN prediction of the primordial abundances of the light elements, the χ lifetime is constrained as τ χ < ∼ 1 sec [23]. Equations (10) and (23) show that the entropy ratio is proportional to the relic density Ω χ h 2 that χ would have if it were stable. We saw in Fig. 2 that for fixed coupling G this quantity is maximal if T F ∼ m χ ; more accurately, the maximum of Ω χ h 2 is achieved for x F = 1.8 (2.1) if χ particles annihilate from an S− (P −)wave initial state. Entropy production by late χ decays is thus most efficient when the χ particles decoupled semi-relativistically, with their lifetime fixed to the maximal value of ∼ 1 sec.
We can construct a feasible scenario that fulfills these conditions by introducing a sterile neutrino which mixes with an ordinary neutrino. Here we treat both m χ and the mixing angle θ as free parameters. In sharp contrast to conventional cosmological scenarios with sterile neutrinos [10], the sterile neutrino is assumed to be in thermal equilibrium in the early universe. In ordinary sterile neutrino models, thermal equilibrium is not reached because the Yukawa coupling of the sterile neutrino with SM particles is tiny. One possible method for the χ pair production and annihilation to reach thermal equilibrium is to extend sterile neutrino models by adding another hypothetical boson Z ′ . Let it have coupling g Z ′ f with the SM fermion pair ff, and g Z ′ χ with the sterile neutrino pair. If the Z ′ -boson mass m Z ′ is larger than the χ energy, the χ annihilation cross section has the form of Eq.(18) Although g Z ′ f and m Z ′ are constrained by high energy experiments, g Z ′ χ can be as large as unity. Therefore, χ annihilation can be in thermal equilibrium before its semi-relativistic decoupling. Decoupling occurred at T ∼ m χ if m χ ≃ 1 GeV · 3 · 10 −9 GeV −2 /G 2/3 g 1/6 * . In order to estimate the amount of entropy released by the decay of sterile neutrinos in this setup, we have to calculate their lifetime. For simplicity we ignore propagator effects. When the sterile neutrino mass is smaller than the W −boson mass m W = 80 GeV, it decays into three SM fermions, with decay width where sin 2 θ W = 0.23 is the weak mixing angle. When the sterile neutrino mass is larger than the Z−boson mass m Z = 91 GeV, the sterile neutrino predominantly decays into a SM gauge boson and a lepton. Its decay width is then proportional to m 3 χ , and given by In the in-between case where m W < m χ < m Z , we obtain Figure 3 shows contours of the entropy increase s f /s i due to sterile neutrino decay in the (1/ √ G, sin θ) plane. We set the freeze-out temperature to x F = 2.1, which maximizes m χ Y χ,i ; this can be achieved by chosing the mass m χ appropriately.
The thick line indicates the BBN limit on the sterile neutrino lifetime, τ χ = 1 sec. Eq. (23) shows that for given neutrino mass, the released entropy will be maximal if θ is chosen such that τ reaches this upper limit.
The behavior of the contours in Fig. 3 is easy to understand from Eq.(23). In the relevant limit θ ≪ 1 and keeping g * constant, we have τ χ ∝ θ −2 m −5 χ (θ −2 m −3 χ ) for m χ < (>) m W . The entropy ratio thus scales as θ −1 m Along the τ χ = 1 sec contour, the entropy release increases proportional to m χ ∝ G −2/3 both for m χ < m W and for m χ > m W . Fig. 3 can be extended to even smaller G, i.e. larger Z ′ masses, so long as m χ is smaller than the re-heat temperature after inflation, so that χ was in thermal equilibrium in the RD epoch. If at the same time θ is decreased so that τ χ = 1 sec remains constant, very large entropy dilution factors could be realized, This result is only valid if the mixing-induced interactions of χ are not in thermal equilibrium for T < ∼ m χ . Since these interactions are also responsible for χ decay, this assumption is satisfied whenever τ χ ≫ t F ; we saw in the discussion of Eq. (23) that this strong inequality is in any case a condition for significant entropy release from χ decay.
We finally note that for given m χ the entropy released in χ decays is maximal if G is so small that χ was ultra-relativistic at decoupling, since this maximizes Y χ,eq (x F ). Again setting τ χ = 1 sec by appropriate choice of θ, this yields s f s i ≤ 10 4 · m χ 10 3 GeV . (28)

Conclusion
In this paper we have developed an approximate analytic method for calculating the thermally-averaged annihilation cross section of semi-relativistically decoupling particles and for estimating their relic density. We have shown that this approximate solution can be smoothly matched to the well-known non-relativistic approximation at the point of intersection. We have argued that such relics cannot form the observed cosmological dark matter. However, we pointed out that the late decay of metastable semi-relativistically decoupling relics can be an efficient source of entropy production. As an example of this entropy production mechanism we discussed a scenario with a sterile neutrino, and illustrated to what extent entropy can be increased.

Acknowledgments
This work was partially supported by the Marie Curie Training Research Network

Appendix A: Modified Bessel Functions
In this appendix, we summarize some properties of the modified Bessel function.
Using an integral representation, the modified Bessel function of the second kind is defined by In particular, the calculation of the relic abundance involves K 1 (z) and K 2 (z), The lower order terms of the series expansion of K 1 (z) and K 2 (z) are given by The asymptotic expansion of K ν (z) is given by Appendix B: Validity of the Maxwell-Boltzmann

Distribution
In the calculations of this paper we used the Maxwell-Boltzmann (MB) distribution also for particles that were semi-relativistic at decoupling; this assumption is e.g.
implicit in Eq. (12). At first sight this seems quite dangerous. For example, at T = m χ , i.e. x = 1, the MB result for n χ,eq overestimates the Fermi-Dirac distribution by about 7%, and underestimates the Bose-Einstein distribution by about 10%. Since χ annihilation always involves two χ particles, one might assume that the total error associated with the use of the MB distribution is about twice as large. In this Appendix we show that the MB distribution can indeed be used to compute the thermally averaged cross section and the decoupling temperature as long as x F > ∼ 1. For smaller x F , one has to use the proper Fermi-Dirac or Bose-Einstein distribution only in the very last step, when calculating Y χ,∞ .
We begin by expanding the true distribution function, where the upper (lower) sign is for fermionic (bosonic) χ particles. Note that the correction term in parentheses has exactly the same form as the "statistics factors" appearing in the collision term of the full Boltzmann equation [2]. For consistency these statistics factors therefore also have to be included. Up to first order in these correction factors, the temperature dependent terms in the integrand defining the collision term for χχ ↔ ff processes then read for fermionic f : here c χ = f χ /f χ,eq is independent of energy as long as χ is in kinetic equilibrium (through elastic scattering on SM particles); in that case we can equivalently write c χ = n χ /n χ,eq . In order to derive the full collision term, I has to be multiplied with the squared matrix element and integrated over phase space [2].
In the usual treatment of WIMP decoupling, all the exponential terms in the square parentheses are neglected, so that the collision term becomes proportional to n 2 χ − n 2 χ,eq times the thermally averaged cross section defined in Eq. (12). Unfortunately the full correction term introduces additional dependence on the final state energies E f and Ef . In order to keep the numerics manageable, we assume that they can be replaced by E χ 1 and E χ 2 , respectively. This is certainly true (by energy conservation) for the sum E f + Ef ; this has already been used in deriving Eq.(34).
Inserting this corrected collision term into the Boltzmann equation, and following the formalism of [7], finally yields a modified thermally averaged cross section times initial state velocity: σv = 1 n 2 χ,eq with F = 2s 1 − 4m 2 χ /s, E + = E χ 1 + E χ 2 and E − = E χ 1 − E χ 2 . This reduces to Eq.(12) if the expression in square parentheses is simply replaced by 1. In the following we assume that χ particles annihilate from an S−wave. P −wave annihilation would favor larger energies, where the correction terms in Eq.(36) are smaller. Note that we also have to use the expanded form (33) of the distribution function when calculating n χ,eq in Eq.(36); otherwise the solution of the Boltzmann equation will not yield n χ ≃ n χ,eq , including the correction terms, at T ≫ T F .   (12), while the dashed (red) curve shows this ratio if Eq.(12) is replaced by our approximation (17). The dotted (blue) curve is the same as the solid curve in Fig. 1. The size of the correction terms in Eq.(36) is illustrated by the solid (black) curve in Fig. 4. We see that the correction amounts to less than 2% for all x > ∼ 1. This is due to a strong cancellation between the corrections in the integrand of Eq.(36) and those in the overall factor 1/n 2 χ,eq . The dashed (red) curve shows that for x ∼ 2 the errors due to the use of the MB distribution and due to our simple parameterization (17) add up, leading to a total error of about 2.7% at most. The Fermi-Dirac corrections to the thermally averaged cross section begin to be significant for x < ∼ 0.5. However, here one enters the ultrarelativistic regime, where the final relic density is no longer sensitive to the decoupling temperature. We therefore expect the effect of using the MB distribution in Eq.(12) on the final prediction of the relic density to be quite small throughout. Ratio Ω χ,MB / Ω χ Ω χ,mix / Ω χ Figure 5: Effect of using the Maxwell-Boltzmann distribution on the predicted relic density, calculated using the approximation Y χ,∞ = Y χ,eq (x F ), for fermionic χ particles annihilating from an S−wave initial state. The dashed (red) curve shows the ratio of the prediction using the Maxwell-Boltzmann distribution everywhere to the corrected prediction based on Eq.(33) and (36). The solid (black) curve shows the analogous ratio, where correct Fermi-Dirac distribution has been used to evaluate Y χ,eq (x F ), but σv and x F have still been obtained using the MB distribution.
Parameters are as in Fig. 2.
This is illustrated in Fig. 5, where the relic density has been calculated using the simple assumption Y χ,∞ = Y χ,eq (x F ); we have used the same parameters as in Fig. 2. The dashed (red) curve shows that using the MB distribution everywhere will overestimate the relic density for m χ < ∼ 5 MeV, i.e. for x F < ∼ 1. However, the black curve shows that this can easily be corrected by using the Fermi-Dirac distribution only in the final step, i.e. when calculating Y χ,eq (x F ); σv and x F can still been calculated using the MB distribution. This validates our treatment in the main text.