Effects of quantum statistics on relic density of dark radiation

The freeze-out of massless particles is investigated. The effects due to quantum statistics, Fermi-Dirac or Bose-Einstein, of all particles relevant for the process are analyzed. Solutions of appropriate Boltzmann equation are compared with those obtained using some popular approximate methods. As an application of general results the relic density of dark radiation in Weinberg’s Higgs portal model is discussed.


Introduction
In recent decades the cosmic microwave background radiation has been extensively analyzed, unveiling many crucial facts about the history of the Universe and its constituents. Although we are convinced that about 27% of the total massenergy fraction of the Universe consists of presumably cold or warm dark matter, there are still some hints that additional form of dark radiation (DR) -called also hot dark matter -may also exist. According to the recent Planck satellite measurements [1] the effective number of light neutrino species N eff varies (depending on the effects included) from 2.99 ± 0.20 to 3.15 ± 0. 23 (1σ ), which indicates that both the Standard Model (SM) and models with fractional N eff ≡ N eff − N SM eff > 0 are consistent with this result (one fully thermalized neutrino i.e. N eff = 1 is ruled out at over 3σ ). On the other hand, the value of the Hubble constant H 0 obtained from direct observations [2] performed using the Hubble Space Telescope i.e. 73.24 ± 1.74 km s −1 Mpc −1 is significantly bigger than 67.8 ± 0.9 km s −1 Mpc −1 inferred from the Planck data. Such high value of H 0 favors additional contribution to N eff (even of order 0.6 [3]), however global fits still prefer the standard CDM scenario [4]. There is hope that this apparent tension will be clarified in near future, especially with the advent of the CMB-S4 experiment [5] that will a e-mail: Marek.Olechowski@fuw.edu.pl b e-mail: Pawel.Szczerbiak@fuw.edu.pl be able to probe N eff with precision better than eq0.03. So precise experimental results will require more accurate theoretical treatment of DR freeze-out details than is usually considered in the literature.
One of the most widely studied scenarios predicting an existence of additional form of radiation is the Weinberg's Higgs portal model [6] in which a massless Nambu-Goldstone boson of a broken global U(1) symmetry gives fractional contribution to N eff . It is often assumed that the freeze-out of that boson takes places just before μ + μ − annihilation, resulting in N eff ≈ 0.39. However, there are some effects, which we discuss in this work, that might strongly influence the Goldstone boson relic density in such case. Moreover, we consider situations when such boson decouples at different temperature. Most of our results may be applied also in other models of DR.
The outline of this work is as follows. In Sect. 2 we discuss the Boltzmann equation describing the freeze-out of a massless particle and some approximate methods used in the literature. Main features of the Weinberg's Higgs portal model are given in Sect. 3. In Sect. 4 the relic density of DR in this model is analyzed in some detail using the general results presented in Sect. 2. Section 5 contains our main conclusions.

Boltzmann equation for massless particles
Boltzmann equation in FLRW metric takes the form where f ( p, t) is a distribution function to be found, whereas terms at the RHS are collision integrals for elastic (E) and inelastic (I ) processes. Elastic collisions are responsible for momentum exchange between particles (preserving their relative numbers) and in consequence help to maintain kinetic equilibrium. Inelastic collisions are related to (co)annihilation processes that influence chemical equilibrium. Elastic processes are usually stronger than inelastic ones, so we will assume here that kinetic equilibrium is maintained. Such assumption is not always valid and under certain circumstances may lead to sizable differences in relic density calculation. There are several methods which are applicable for Boltzmann equation solution without kinetic equilibrium -this problem has been extensively analyzed in the context of neutrino decoupling in the Standard Model. We can distinguish two major classes of relic density calculation for relativistic particles in such a case: discretization in momentum space [7][8][9][10] and spectral methods based on fixed [11,12] or dynamical [13,14] basis of orthogonal polynomials. All these approaches allow for high accuracy but at the price of lengthy and complicated expressions that are unpractical for qualitative study. In the present work we are mainly interested in estimation of Bose-Einstein (BE)/Fermi-Dirac (FD) statistics corrections to the relic density for massless particles. Thus, it suffices in our case to perform analysis assuming the kinetic equilibrium.
The collision integral for inelastic processes involving particle χ , C I ( p χ , t), with only two-body processes χa ↔ bc taken into account, may be expressed as where the sum over {a, b, c} corresponds to all allowed processes χa ↔ bc. We will calculate the above collision integral exploiting the method proposed in [15], where the effect of FD statistics was analyzed in the context of relic density calculation for relativistic fermions (0 = m T f ) staying in kinetic equilibrium. Below we will focus on annihilation of massless particles (bosons and fermions) including full statistics for both initial and final states. The results will be used in Sect. 4 to analyze the Goldstone bosons freeze-out in the Weinberg's Higgs portal model.
Let us consider annihilation process of the formχχ → N N, where N stays in kinetic and chemical equilibrium during χ freeze-out. For definiteness, we will focus on one process of this kind (generalization is straightforward). Distribution functions for χ and N may be written in the form where z = z(t) is the so-called chemical pseudopotential (equal for χ andχ ) [16]. We also defined x ≡ m/T and y ≡ E/m (for one massive annihilation product N it is convenient to choose m ≡ m N = 0). The statistical factors s in and s out equal +1 and −1 for fermions and bosons, respectively -by setting these factors to zero one obtains the Maxwell-Boltzmann (MB) approximation. Using the above definitions and integrating expression (2) over χ momenta one may rewrite Eq. (1) in the following form [15] dz dx where g is the number of χ degrees of freedom, the Hubble parameter is given by and 1 J n (z, x) ≡ ∞ 0 dy y n e xy (e xy+z + s in ) 2 .
The most elaborate part is the calculation of the collision integral S I which, after using four-momentum conservation, is 5-dimensional. We use the following convenient set of variables: energies of the incoming particles, E 1 and E 2 , one of the outgoing particles' energy, E 3 , the angle between momenta of the incoming particles, θ , and the acoplanarity angle φ. Defining dimensionless parameters u ≡ E 1 /m, v ≡ E 2 /m and t ≡ E 3 /m one can write where dt × 1 e xu+z + s in e xv+z + s in (e xt + s out ) e x(u+v−t) + s out (10) and we introduced the following functions of the parameters: If (which is often the case in leading approximation) M χχ→NN depends only on s (given here by s = 2uv m 2 (1 − cos θ)), introducing new combinations of the parameters p = x(u + v) and q = x √ u 2 + v 2 + 2uv cos θ , we can reduce the above integral to In such a case the integration over φ in Eq. (9) is trivial and the expression simplifies. Then, the Boltzmann equation (5) may be written in the following integro-differential form Note that we introduced dimensionless (in contrast to S I (z, x)) collision integral (see Eq. (11)) The number density of χ (=χ ) may be easily found after integration over the whole energy range 2 2.1 MB approximation and limited inclusion of BE/FD statistics Let us consider the particle number density normalized by the entropy density: Y (x) ≡ n(x) s(x) . In chemical and kinetic 2 If χ andχ are distinguishable one should consider n = n χ + nχ instead. equilibrium (subscript eq) we just have (using Eq. (14) with where For the Maxwell-Boltzmann statistics one can write Then, the Boltzmann equation (12) simplifies to (we use s in = s out = 0 in Eqs. (8) and (10)) One can also rewrite it in the Lee-Weinberg form where v is the Möller velocity for the incoming particles and Please note the presence of additional factor 1/ζ in in the RHS of Eq. (19) as compared to the standard form [17]. It comes from the fact that Y eq contains information about the incoming particles statistics whereas Eq. (12), after putting everywhere s out = s in = 0, does not. MultiplyingS I (z, x) by ζ in one can obtain the familiar Lee-Weinberg equation. Thus, we have two options: 1. Use Eqs. (15) and (19) with ζ in = 1 (pure MB approximation). 2. Use Eqs. (15) and (19) with ζ in given by Eq. (16) in the former but ζ in = 1 in the latter (we call it fractional inclusion of BE/FD statistics for incoming particles and denote as fBE/fFD).
One can also include BE/FD statistics (again, only for incoming particles) in thermally averaged cross section where fh is a hyperbolic function depending on the statistics of incoming particles: fh ≡ sinh (cosh) for bosons (fermions) and index p stands for partial inclusion of statistics. Thus, we consider another approximate method: 3. The same rules as in the point 2. above but with σ v p (Eq. (21)) instead of σ v MB (Eq. (20)) (we call it partial inclusion of BE/FD statistics for incoming particles and denote as pBE/pFD).
The main difference between Eqs. (11) and (21) is that in the latter case the expression does not depend on Y (or equivalently on z), which significantly simplifies calculations. Note also that none of the above approximations include any effect of the outgoing particles statistics.

Weinberg's Higgs portal model
For easy reference we remind here the main features of the model proposed in [6] (using notation mainly from [18]). In addition to the SM fields it contains a dark sector consisting of a complex scalar φ and a Dirac fermion ψ. Both new fields are charged only under a global U (1) dark symmetry with charges: A non-zero VEV of φ spontaneously breaks this symmetry to the discrete Z 2 parity. The Goldstone boson associated with this symmetry breaking contributes to DR while the lighter dark fermion (ψ gives two Majorana fermions with different masses) is stable and plays the role of cold dark matter (CDM). We will concentrate here on the scalar sector which is the most important one for our analysis. The corresponding part of the Lagrangian density is given by where H is the Standard Model Higgs doublet. After the spontaneous symmetry breaking the scalar fields may be written as where σ is the Goldstone boson of the broken U (1) dark and contributes to DR. Two neutral scalarsh andρ mix, forming two mass eigenstates (h and ρ): with masses and the mixing angle given by the condition The scalar potential in (22) depends on five parameters. They must satisfy two conditions in order to give the correct values of v H and the mass of the Higgs particle (identified with h).
As the three parameters, describing the remaining freedom in the scalar Lagrangian (22), we choose two couplings, λ φ and κ, and the mass of the non-SM-like scalar, m ρ . Other parameters present in the mass formulas (25) and (26) may be expressed as where the approximate equalities are valid if the scalar ρ is relatively light and does not mix strongly with the SM-like Higgs i.e. m 2 ρ m 2 h , sin 2 θ 1. Let us now estimate phenomenologically interesting range of values for κ and λ φ that will be used in our numerical scan. The contribution from the dark particles to the invisible Higgs decay width may be approximated as where we conservatively neglected the Higgs decays into dark fermions. The LHC constraint on the invisible Higgs decays i.e.
where B inv 24% [19] and h,SM 4 MeV, may be used together with Eq. (31) to obtain the following upper bound on the portal coupling: Bounds on the singlet scalar self-coupling, λ φ , are much weaker. They follow from the requirement of perturbativity to some high energy scale and may be obtained from the RG equations [20]: where t ≡ Q 2 /Q 2 0 . The scalar couplings stay perturbative up to the GUT energy scale if at low energy λ φ 0.25. On the other hand, if λ φ = 1 at the weak scale, it becomes nonperturbative already at scale of order 100 TeV. Both values i.e. λ φ = 0.25, 1 will be our reference points in numerical scans in the next section.

Dark radiation in Weinberg's Higgs portal model
In this section we will compare solutions of the Boltzmann equation (12) with the instantaneous freeze-out approximation as well as with three approximate methods introduced in Sect. 2.1. The model of DR proposed by Steven Weinberg [6] will be used as a testing ground for this purpose.
In order to calculate the relic abundance of the Goldstone boson σ one needs to know its cross section for the annihilation into SM particles. In the original paper [6] it was assumed that σ decouples at temperature (just) above the muon mass. Then, in some analyses (e.g. [18,21]) only the annihilation of σ into muons was taken into account. However, pions are not much heavier than muons and their contribution should also be included. It was shown [22,23] that the effect from pions may be a few times stronger than the one from muons (in the context of Weinberg's Higgs portal model it was analyzed e.g. in [24]). The squares of matrix elements for the relevant processes with ρ exchanged in the s channel may be approximated as 3 We checked that the effects from τ ± , gluons and free quarks become important at temperatures T 0.4 ÷ 0.5 GeV. This suggests that the Goldstone bosons stay in equilibrium at temperatures above QCD . Thus, we assume in our analysis that σ particles indeed are in equilibrium for T > 0.4 GeV. It occurs that from the phenomenological point of view the most interesting region of the parameter space is that close to the resonance i.e. s ∼ m 2 ρ . If the resonance is narrow, i.e. when ρ /m ρ 2 1, one may use the approximation We assume for simplicity that the CDM fermion is not very light and fulfills the condition m CDM > m ρ /2. Then the width of ρ is dominated by its decays into Goldstone bosons pairs, giving Because ρ /m ρ 2 = (λ φ /16π 2 ) 2 1 approximation (37) may be used even for λ φ as big as 1. We do not consider larger values of λ φ because, as explained below Eq. (34), they would lead to a cut-off scale below 100 TeV. Thus, we may apply the narrow-resonance approximation (37) and rewrite the matrix elements (35) and (36) in the following form Using this approximation we can also rewrite integrals (11), (20) and (21) as σ v MB = spins |M| 2 512π which considerably simplifies numerical calculations. 4 Before presenting the results of our numerical analysis let us remind the definition of the effective number of neutrino species, N eff . It is usually defined by the formula where ρ R , ρ γ and ρ ν correspond to the present total, photons and neutrinos radiation energy density, respectively. Index ifo denotes instantaneous freeze-out approximation. In general we can write where ρ X is an additional form of radiation existing in certain extensions of the Standard Model. In the case of Weinberg's Higgs portal model: X = σ . N eff is normalized in such a way that using the instantaneous freeze-out approximation for neutrino decoupling and putting ρ X = 0 we get N eff = 3. In the Standard Model due to the lack of perfect kinetic 4 Equation (41) follows from (11) after changing variables and performing one integration, using (37). The functions appearing in the integral (41) depend on the integration variable y as follows: p = x(m ρ /m μ ± )y, q = x(m ρ /m μ ± ) y 2 − 1 and V = 1 − 4(m μ ± /m ρ ) 2 .
equilibrium during neutrino decoupling one gets N SM eff ≈ 3.045 [25]. For simplicity, we will define additional effective number of neutrino species as which for the Goldstone boson σ (additional factor 7 8 due to different statistics and number of degrees of freedom as compared to ν) gives where x end ≈ 20 corresponds to T end ≈ 5 MeV (temperature before neutrino decoupling when the freeze-out process of σ is completed).

Instantaneous freeze-out approximation
Instantaneous freeze-out approximation is one of the simplest methods for relic density estimation. By definition, freezeout takes place at such x f for which the freeze-out parameter η(x), being the ratio between the interaction rate and the Hubble parameter H , equals one: At temperatures below T f (x > x f ) the yield of considered particles remains constant i.e. Y (x) = Y eq (x f ). For σ v = σ v MB (see Eq. (20)) with matrix elements squared as in (39)-(40) and ρ given by (38), one gets the following condition for where In our numerical calculations we use the number of the effective degrees of freedom, g * (x), as given in [27]. In the upper panels of Fig. 1 we plotted N eff obtained in the instantaneous freeze-out approximation, Y = Y eq (x f ), as a function of m ρ and κ (blue lines), using definition (47) and the freeze-out temperature calculated from Eq. (49). The results are compared to those obtained by solving the The freeze-out process for parameters corresponding to the black dots for κ = 2 · 10 −4 is presented in more detail in Fig. 2 Boltzmann equation (red lines) -see Sect. 4.2. In almost all cases the instantaneous freeze-out approximation overestimates N eff for a given values of κ and λ φ (dark blue lines are below corresponding red lines). The best agreement between the instantaneous freeze-out approximation and the full result is achieved for N eff close to 0.1, although in that region the instantaneous freeze-out approximation works only for m ρ 1 GeV. Note also that neglecting the σ annihilation into pions (given by Eq. (40)), while keeping only annihilation into muons (39) (light blue lines in Fig. 1), leads to quite substantial underestimation of the resulting N eff (the difference amplifies when κ decreases). 5 It is worth emphasizing that the instantaneous freeze-out approximation breaks down for very small κ and m ρ 1 GeV. Thus, exploration of this part of the parameter space requires a more accurate approach, which we will discuss in the next section. Let us only add that the above-mentioned problem is related to the stiffness of condition (48). As we can see in the lower panels of Fig. 1, for small enough κ the η(T ) function never reaches 1.
The current experimental limits, shown in Fig. 1, exclude κ 10 −2 and the region of small m ρ with κ 10 −3 . Note also that we conservatively do not consider here the CDM sector, which could (depending on the dark matter mass) constrain the parameter space even further [21,26].

Solutions of the Boltzmann equation
In this subsection we compare the results obtained from solutions of the Boltzmann equation (with statistics of incoming and outgoing particles taken into account) with those obtained with the approximate methods introduced in Sect. 2.1 (which include some effects of incoming particles statistics) and with the instantaneous freeze-out approximation discussed in the previous subsection.
Regions of small κ, where the last approximation breaks down, are especially important from the viewpoint of near future experiments that will be able to probe N eff with accuracy even better than 0.03. Let us start the discussion by considering two sample points from the parameter space with small κ: two lower (black) dots in the upper left panel of Fig. 1.
In Fig. 2 we show corresponding evolutions of Y (T ) and z(T ) for all three approximations described in Sect. 2.1 as well as for the most accurate of our methods. Line colors i.e. blue, violet, green and red correspond to the pure MB approximation (point 1. under Eq. (20)), fractional inclusion of incoming particles statistics (point 2. under Eq. (20)), partial inclusion of incoming particles statistics (point 3. under Eq. (21)) and full inclusion of particles statistics (Eq. (12)), respectively. As we have already seen in Fig. 1, the instantaneous freeze-out approximation (dashed blue line in the upper right panel in Fig. 2) overestimates the relic density, as compared to the solutions of the Boltzmann equation (with full or any approximate inclusion of statistics effects). Such behavior is related to the convexity of g * (T ) function, which for T ∼ 100 ÷ 200 MeV and 30 ÷ 50 MeV is characterized by strong variability. As a result, it is harder for σ particles to follow the equilibrium density i.e. larger κ is required as compared to the instantaneous freeze-out approximation, where the effect from g * (T ) is point-wise. We checked that the instantaneous freeze-out approximation gives results closer to those obtained by integrating the Boltzmann equation when g * (T ) during the σ freeze-out changes relatively mildly (e.g. for T 200 MeV). However, the accuracy obtained with the instantaneous freeze-out approximation is typically worse than the accuracy of other approximate methods discussed in this work. One can also include the backreaction of g σ (x) on g * (x), however we checked numerically that this effect is negligible in the whole range of the analyzed parameter space.
It is also worth noting in the left panels of Fig. 2 that σ freezes in for T ∼ 60 ÷ 130 MeV. For such temperatures the pseudopotential z(T ) decreases with time, which is related to the fact that the annihilation cross section for small m ρ reaches its maximum in this temperature range (see also η(T ) dependence in the lower panels of Fig. 1).
Having discussed general features of the lines presented in Fig. 2, let us now take a closer look at the differences between approximations in solving the Boltzmann equation, described in Sect. 2.1. We checked both analytically and numerically that in the range of analyzed parameter space the phase space integral (41) may be expressed for different configurations of incoming/outgoing particles statistics as 6 where (see Eq. (42)) This clearly shows that for the cases when the incoming and outgoing particles have different statistics (BEFD or FDBE) the effects of statistics of initial and final states cancel each other to large extend and the MB approximation works quite well. When the initial and final particles have the same statistics one gets amplification (BEBE) or suppression (FDFD) with respect to the MB approximation. In Weinberg's Higgs portal model one has to consider a combined case i.e. BEFD (annihilation into muons) and BEBE (annihilation into pions). The effect from the latter (factor coth 2 m ρ x/4m μ ± ) starts to be important for small m ρ and x (larger T ) and also for small κ. The dependence on κ follows from the fact that σ with smaller κ decouples at higher temperature i.e. at smaller x. These effects can be seen in Fig. 3, where for small enough m ρ and κ red curves are placed partially under the blue ones i.e. in order to get a given value of N eff smaller κ is needed when using the Boltzmann equation with full statistics effects as compared to the MB approximation. One can also observe from Fig. 3 that for a given value of κ and m ρ in most cases N pBE eff < N MB eff < N fBE eff and relative differences between approximations decrease when κ increases. Only for small enough κ and m ρ we have N pBE eff > N MB eff . These relations can be easily understood with the help of Eqs. (51) and (52). Moreover, for the most part of the parameter space (i.e. except the region with sufficiently small m ρ and κ, where pions statistics matters) full inclusion of incoming and outgoing statistics results in smaller N eff than the above-mentioned approximations. There are two effects which explain this behavior. Firstly, from Eq. (14) one can see that for given z and x inclusion of incoming particles statistics (here: BE) results in n(x) (equivalently Y (x)) smaller than in the MB approximation (see Eq. (17)) and the difference increases with z. Secondly, comparing Eqs. (12) and (18) one can show (see Definition (8)) that for given z and x coefficients multiplying two terms in the RHS of Eq. (12) (J 3 (z, x)/J 2 (z, x) and 1/J 2 (z, x) respectively) are smaller than those in Eq. (18). First term in the RHS of Eq. (12) dominates for small z (when Y (x) traces Y eq (x)) which causes weaker change in |dz/dx| as compared to the MB case. When the second term starts to be important (for larger z) the coefficient 1/J 2 (z, x) effectively weakens the annihilation strength (given byS I (z, x)) leading to faster |dz/dx| change with respect to other approximations. Both effects are visible in Fig. 2. Contrary to naive expectations the MB approximation gives better accuracy than the fBE one, which includes the effect from incoming particles statistics only in Y eq . In order to obtain more precise results, it is necessary to take into account the effects from statistics also in the calculation of the thermal average of the appropriate cross section (pBE). Let us add that taking into account σ annihilation into muons only (see Figs. 1, 3) leads to sizable discrepancies with respect to the case when full annihilation is considered. The only exception holds for m ρ ≤ 2m π (i.e. when annihilation into pions does not take place) -one can observe that continuous lines rapidly change slope and red lines converge towards dashed ones.
Relative differences between our approximations may reach N eff ∼ 0.05 and more. 7 Thus, statistics of both incoming and outgoing particles are relevant, especially for moderate values of κ and m ρ . Minimal value of N eff obtained in the scan (∼ 0.06) is related to the starting moment for the calculation i.e. T start = 400 MeV, however for larger T Y eq (T ) does not change significantly. Therefore, if σ freezes out during or after the QCD phase transition, it shall give contribution to N eff measurable by near future experiments. On the other hand, maximal N eff equals approximately 0.5 and is achieved for large κ and m ρ near the h 1 resonance (see Eq. (37)), which is not preferred due to collider and astrophysical bounds. It is worth mentioning that the widely studied scenario with N eff ∼ 0.39 for moderate values of λ φ can be probed in major parts of the parameter space by the ILC (and new generation of CMB satellite experiments).
However, the region with smaller λ φ may be more challenging in this context as the lines of constant N eff move towards smaller κ, becoming harder to be probed.
Let us finish this section with the discussion of differences between possible combinations of incoming and outgoing particles statistics (see Eqs. (53) and (54)). In Fig. 4 we show similar plots to those presented in Fig. 3 but for 8 m ρ ≤ 1.2 GeV, κ ≤ 10 −3 and taking into account only (dominant) annihilation into pions (Eq. (40)).
Then, BEBE case corresponds to Weinberg's Higgs portal model with annihilation into muons neglected, while the others correspond to toy models with statistics of incoming and/or outgoing particles changed to the FD one. For mixed statistics BEFD/FDBE continuous red curves are placed above/below those that were obtained including the effects from incoming particles statistics only i.e. fBE/fBE and pBE/pFD (as well as from the MB approximation). For such mixed cases the effects of statistics of incoming and outgoing particles approximately cancel out in D (see Eq. (53)). The main reason for smaller/bigger values of N eff in such mixed cases, especially when compared to the results of the MB approximation, is the presence of the statistics factor s in = ± 1 in Eqs. (8) and (14) (as discussed below Eq. (55)). As one can see (and what can be observed also in Fig. 4) for homogeneous statistics (BEBE and FDFD) the additional factors (coth 2 m ρ x/4m μ ± and tanh 2 m ρ x/4m μ ± ) in D given by (54)) are important for small exchanged particle mass (here m ρ ) and coupling (here κ). The effects of outgoing particle statistics decrease with m ρ and κ and eventually become negligible.

Conclusions
We have investigated the problem of calculating the relic abundance of Dark Radiation which freezes out before the SM neutrinos decoupling. We used the Boltzmann equation for relativistic particles with their statistics taken into account. This method was compared to the instantaneous freeze-out approximation and to some approximate methods in which statistics of DR particles is included only in a limited way or completely ignored. As an interesting illustration of all these methods we analyzed in some detail the relic density of DR -measured by the change of the effective number of neutrino species N eff -in the Weinberg's Higgs portal model. The main results are as follows: • The popular instantaneous freeze-out approximation can not be applied for small values of m ρ and κ for which the Boltzmann equation must be used. • In most of the remaining regions of the parameter space the instantaneous freeze-out approximation overestimates N eff . The main reason is convexity of g * (T ) which in this simple method is not taken into account.
The instantaneous freeze-out approximation gives best results for N eff ∼ 0.1 and m ρ 1 GeV. The resonant exchange of the scalar ρ (for m ρ up to a few hundreds MeV) also plays an important role. • When calculating the annihilation cross section of DR particles (Goldstone bosons σ ) it is crucial to include muons and pions as the final states. Pions are quite often ignored which may lead to underestimation of N eff by as much as 0.1.
• Not taking (fully) into account the statistics of DR particles and particles into which DR annihilates may change the obtained values of N eff by up to about 0.05. Contrary to naive expectation, in some parts of the parameter space ignoring the effects of statistics (MB approximation) may give better prediction for N eff than inclusion of only some of the effects -those in evaluation of Y eq (T ) -due to the DR statistics (fBE approximation). Inclusion of more effects from the DR statistics (pBE approximation) leads to more accurate results than those obtained using the simplest MB approximation.
The present experimental data leave quite substantial uncertainty in determining the value of N eff . One may expect that results of near future experiments will lead to much better determination of N eff and will allow to test scenario considered in this work. In such a case it will be important not only to use the statistics of incoming and outgoing particles in the appropriate Boltzmann equation but maybe also to take into account the effects of deviations from kinetic equilibrium. This may be especially important for cases with resonant exchange of ρ when DR particles σ may decouple kinetically before thermal decoupling. This might happen because in elastic scattering particle ρ is exchanged in the t channel for which there is no resonant enhancement [28]. Another potentially important issue is the relation between the DR and CDM sectors. We assumed in our analysis that CDM particles are so heavy that they do not influence the DR properties. However, it may occur, especially when future more precise experimental results are available, that the case of lighter CDM should be considered simultaneously with DR.