Effects of quantum statistics on relic density of dark radiation

The European Physical Journal C, Sep 2018

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.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

https://link.springer.com/content/pdf/10.1140%2Fepjc%2Fs10052-018-6195-0.pdf

Effects of quantum statistics on relic density of dark radiation

The European Physical Journal C September 2018, 78:704 | Cite as Effects of quantum statistics on relic density of dark radiation AuthorsAuthors and affiliations Marek OlechowskiPaweł  Szczerbiak Open Access Regular Article - Theoretical Physics First Online: 01 September 2018 Received: 23 July 2018 Accepted: 27 August 2018 52 Downloads Abstract 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. 1 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 mass-energy 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_{\mathrm{eff}}\) varies (depending on the effects included) from \(2.99\pm 0.20\) to \(3.15\pm 0.23\;(1\sigma )\), which indicates that both the Standard Model (SM) and models with fractional \(\Delta N_{\mathrm{eff}}\equiv N_{\mathrm{eff}}-N_{\mathrm{eff}}^\mathrm{SM}>0\) are consistent with this result (one fully thermalized neutrino i.e. \(\Delta N_{\mathrm{eff}}=1\) is ruled out at over 3\(\sigma \)). 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\pm 1.74\;\mathrm km\,s^{-1}Mpc^{-1}\) is significantly bigger than \(67.8\pm 0.9\;\mathrm km\,s^{-1}Mpc^{-1}\) inferred from the Planck data. Such high value of \(H_0\) favors additional contribution to \(\Delta N_{\mathrm{eff}}\) (even of order 0.6 [3]), however global fits still prefer the standard \(\Lambda \)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 be able to probe \(\Delta N_{\mathrm{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_{\mathrm{eff}}\). It is often assumed that the freeze-out of that boson takes places just before \(\mu ^+\mu ^-\) annihilation, resulting in \(\Delta N_{\mathrm{eff}}\approx 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. 2 Boltzmann equation for massless particles Boltzmann equation in FLRW metric takes the form $$\begin{aligned} E(\partial _t-pH\partial _p)f(p,t)=C_E(p,t)+C_I(p,t), \end{aligned}$$ (1) 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 \(\chi \), \(C_I(p_\chi ,t)\), with only two-body processes \(\chi a\leftrightarrow bc\) taken into account, may be expressed as Open image in new window (2) where the sum over \(\{a,b,c\}\) corresponds to all allowed processes \(\chi a\leftrightarrow 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\not =m\ll 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 \({\bar{\chi }}\chi \rightarrow {\bar{N}}N\), where N stays in kinetic and chemical equilibrium during \(\chi \) freeze-out. For definiteness, we will focus on one process of this kind (generalization is straightforward). Distribution functions for \(\chi \) and N may be written in the form $$\begin{aligned} f_\chi (p,t)\simeq & {} \left( e^{E/T+z}\pm 1\right) ^{-1}\equiv \left( e^{xy+z} + s_{\mathrm{in}}\right) ^{-1}, \end{aligned}$$ (3) $$\begin{aligned} f_N(p,t)\simeq & {} \left( e^{E_N/T} + s_{\mathrm{out}}\right) ^{-1}, \end{aligned}$$ (4) where \(z=z(t)\) is the so-called chemical pseudopotential (equal for \(\chi \) and \({{\bar{\chi }}}\)) [16]. We also defined \(x\equiv m/T\) and \(y\equiv E/m\) (for one massive annihilation product N it is convenient to choose \(m\equiv m_N\not =0\)). The statistical factors \(s_{\mathrm{in}}\) and \(s_{\mathrm{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 \(\chi \) momenta one may rewrite Eq. (1) in the following form [15] $$\begin{aligned} \frac{\mathrm{d}z}{\mathrm{d}x}\simeq \left( A(z,x)\frac{x}{T}\frac{\mathrm{d}T}{\mathrm{d}t}\right) ^{-1}\left( S_I(z,x)-B(z,x)\right) , \end{aligned}$$ (5) where $$\begin{aligned} A(z,x)&\equiv \frac{g}{2\pi ^2}m^3e^z J_2(z,x),\nonumber \\ B(z,x)&\equiv \frac{g}{2\pi ^2}m^3e^zx J_3(z,x)H(x) \left( 1-\frac{1}{1-\frac{x}{3}\frac{g_{*s}'(x)}{g_{*s}(x)}}\right) ,\nonumber \\ S_I(z,x)&\equiv \int \frac{\mathrm{d}^3p}{(2\pi )^3}\frac{1}{E}C_I(p,t), \end{aligned}$$ (6) g is the number of \(\chi \) degrees of freedom, the Hubble parameter is given by $$\begin{aligned} H(x)= \sqrt{\frac{4\pi ^3}{45}}\sqrt{g_{*}(x)}\frac{m^2}{M_{\mathrm{Pl}}}\frac{1}{x^2}, \end{aligned}$$ (7) and1 $$\begin{aligned} J_n(z,x)\equiv \int _0^\infty \mathrm{d}y\,y^n\frac{e^{xy}}{(e^{xy+z}+ s_{\mathrm{in}})^2}. \end{aligned}$$ (8) 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, \(\theta \), and the acoplanarity angle \(\phi \). Defining dimensionless parameters \(u\equiv E_1/m\), \(v\equiv E_2/m\) and \(t\equiv E_3/m\) one can write $$\begin{aligned} S_I(z,x)= & {} \frac{m^4}{512\pi ^6}\left( e^{2z}-1\right) \int {\mathscr {D}}\Phi \int _0^{2\pi }\mathrm{d}\phi \nonumber \\&\times \sum _{\mathrm{spins}}|M_{\chi {\bar{\chi }}\rightarrow N{\bar{N}}}(u,v,t,\cos \theta ,\phi )|^2, \end{aligned}$$ (9) where $$\begin{aligned} {\mathscr {D}}\Phi&\equiv \int _0^\infty \mathrm{d}u\int _0^\infty \mathrm{d}v\int _{-1}^1\mathrm{d}\cos \theta \frac{uv\,e^{x(u+v)}}{\kappa (u,v,\theta )}\int _{t_-}^{t_+}\mathrm{d}t\nonumber \\&\quad \times \frac{1}{ \left( e^{xu+z}+ s_{\mathrm{in}}\right) \left( e^{xv+z}+ s_{\mathrm{in}}\right) \left( e^{xt}+ s_{\mathrm{out}}\right) \left( e^{x(u+v-t)}+ s_{\mathrm{out}}\right) } \end{aligned}$$ (10) and we introduced the following functions of the parameters: \(\kappa (u,v,\theta )\equiv (u^2+v^2+2\,uv\cos \theta )^{1/2}\), \(V(u,v,\theta )\equiv (1-4m^2/s)^{1/2}\), \(t_\pm \equiv \frac{1}{2}\left[ u+v\pm \kappa (u,v,\theta )V(u,v,\theta )\right] \). If (which is often the case in leading approximation) \(M_{\chi {\bar{\chi }}\rightarrow N{\bar{N}}}\) depends only on s (given here by \(s=2uv\,m^2(1-\cos \theta )\)), introducing new combinations of the parameters \(p=x(u+v)\) and \(q=x\sqrt{u^2+v^2+2uv\cos \theta }\), we can reduce the above integral to $$\begin{aligned} {\mathscr {D}}\Phi&= \frac{1}{x^4} \int _{2x}^\infty \mathrm{d}p\int _0^{\sqrt{p^2-4x^2}} \mathrm{d}q\; \frac{1}{(1-e^{-p})(e^{p+2z} - 1)}\nonumber \\&\quad \times \ln \left[ \frac{\cosh \left( \frac{1}{2}(p+q)+z\right) + s_{\mathrm{in}}}{\cosh \left( \frac{1}{2}(p-q)+z\right) + s_{\mathrm{in}}} \right] \nonumber \\&\quad \times \ln \left[ \frac{\cosh \left( \frac{1}{2}\big (p+qV(p,q)\big )\right) + s_{\mathrm{out}}}{\cosh \left( \frac{1}{2}\big (p-qV(p,q)\big )\right) + s_{\mathrm{out}}} \right] , \end{aligned}$$ (11) where now \(V=\left( 1-\frac{4x^2}{p^2-q^2}\right) ^{1/2}\) and \(s=\frac{p^2-q^2}{x^2}m^2\). In such a case the integration over \(\phi \) in Eq. (9) is trivial and the expression simplifies. Then, the Boltzmann equation (5) may be written in the following integro-differential form $$\begin{aligned} \frac{\mathrm{d}z}{\mathrm{d}x}= & {} -\frac{x}{J_2(z,x)} \left( \frac{1}{3}\frac{g'_{*s}(x)}{g_{*s}(x)}J_3(z,x)+ \frac{\sqrt{45}}{256\pi ^{11/2}}\frac{M_{\mathrm{Pl}}}{m} \right. \nonumber \\&\left. \times \frac{\sinh z}{g\sqrt{g_{*}(x)}} \left( 1-\frac{x}{3}\frac{g'_{*s}(x)}{g_{*s}(x)}\right) {\tilde{S}}_I(z,x) \right) . \end{aligned}$$ (12) Note that we introduced dimensionless (in contrast to \(S_I(z,x)\)) collision integral (see Eq. (11)) $$\begin{aligned} {\tilde{S}}_I(z,x)= 2\pi \int {\mathscr {D}}\Phi \sum _{\mathrm{spins}}|M_{\chi {\bar{\chi }}\rightarrow N{\bar{N}}}(s)|^2. \end{aligned}$$ (13) The number density of \(\chi \) (\(={\bar{\chi }}\)) may be easily found after integration over the whole energy range2 $$\begin{aligned} n(x)\equiv n_\chi (x)=g\frac{m^3}{2\pi ^2}\int _0^\infty \frac{y^2\mathrm{d}y}{e^{xy+z}+ s_{\mathrm{in}}}. \end{aligned}$$ (14) 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)\equiv \frac{n(x)}{s(x)}\). In chemical and kinetic equilibrium (subscript eq) we just have (using Eq. (14) with \(z=0\)) $$\begin{aligned} Y_{\mathrm{eq}}(x)=g\,\frac{45}{2\pi ^4}\frac{\zeta ^{\mathrm{in}}}{g_{*s}(x)}, \end{aligned}$$ (15) where $$\begin{aligned} \zeta ^{\mathrm{in}}\equiv \zeta (3) {\left\{ \begin{array}{ll} 1 &{} \text {BE}\\ 3/4 &{} \text {FD} \end{array}\right. }. \end{aligned}$$ (16) For the Maxwell-Boltzmann statistics one can write $$\begin{aligned} Y(x)=e^{-z}Y_{\mathrm{eq}}(x). \end{aligned}$$ (17) Then, the Boltzmann equation (12) simplifies to (we use \(s_{\mathrm{in}}=s_{\mathrm{out}}=0\) in Eqs. (8) and (10)) $$\begin{aligned} \frac{\mathrm{d}z}{\mathrm{d}x}= & {} \ -\frac{g'_{*s}(x)}{g_{*s}(x)}- \frac{\sqrt{45}}{\pi ^{7/2}}\frac{M_{\mathrm{Pl}}m}{x^2}\frac{\sinh z}{g\sqrt{g_{*}(x)}}\nonumber \\&\times \left( 1-\frac{x}{3}\frac{g_{*s}'(x)}{g_{*s}({\tilde{x}})}\right) \langle \sigma v\rangle _{\mathrm{MB}}. \end{aligned}$$ (18) One can also rewrite it in the Lee-Weinberg form $$\begin{aligned} Y'(x)= & {} -\sqrt{\frac{\pi }{45}} \frac{g_{*s}(x)}{\sqrt{g_{*}(x)}} \frac{M_{\mathrm{Pl}}m}{x^2}\nonumber \\&\times \left( Y^2(x)-Y^2_{\mathrm{eq}}(x)\right) \langle \sigma v\rangle _{\mathrm{MB}}\frac{1}{\zeta ^{\mathrm{in}}}, \end{aligned}$$ (19) where v is the Möller velocity for the incoming particles and $$\begin{aligned} \langle \sigma v\rangle _{\mathrm{MB}}= & {} \frac{1}{512\pi }\frac{x^5}{m^5} \int _{4m^2}^\infty \mathrm{d}s\,\sqrt{s}\nonumber \\&\times \sqrt{1-\frac{4m^2}{s}}\,\sum _{\mathrm{spins}}|M(s)|^2\mathrm{K}_1 \left( \frac{x\sqrt{s}}{m}\right) . \end{aligned}$$ (20) Please note the presence of additional factor \(1/\zeta ^{\mathrm{in}}\) in the RHS of Eq. (19) as compared to the standard form [17]. It comes from the fact that \(Y_{\mathrm{eq}}\) contains information about the incoming particles statistics whereas Eq. (12), after putting everywhere \(s_{\mathrm{out}}=s_{\mathrm{in}}=0\), does not. Multiplying \({\tilde{S}}_I(z,x)\) by \(\zeta ^{\mathrm{in}}\) one can obtain the familiar Lee-Weinberg equation. Thus, we have two options: 1. Use Eqs. (15) and (19) with \(\zeta ^{\mathrm{in}}=1\) (pure MB approximation).   2. Use Eqs. (15) and (19) with \(\zeta ^{\mathrm{in}}\) given by Eq. (16) in the former but \(\zeta ^{\mathrm{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 $$\begin{aligned} \langle \sigma v\rangle _{\mathrm{p}}&= \frac{1}{512\pi (\zeta ^{\mathrm{in}})^2}\frac{x^5}{m^5} \int _{4m^2}^\infty \mathrm{d}s\,\sqrt{1-\frac{4m^2}{s}}\,\sum _{\mathrm{spins}}|M(s)|^2\nonumber \\&\quad \times \int _{\sqrt{s}}^\infty \mathrm{d}E_+\frac{e^{-\frac{x}{2m}E_+}}{\sinh \left( \frac{x}{2m}E_+\right) }\nonumber \\&\quad \times \ln \left[ \frac{\mathrm{fh}\left( \frac{x}{4m}\left( E_++\sqrt{E_+^2-s}\right) \right) }{\mathrm{fh}\left( \frac{x}{4m}\left( E_+-\sqrt{E_+^2-s}\right) \right) } \right] , \end{aligned}$$ (21) where \(\mathrm{fh}\) is a hyperbolic function depending on the statistics of incoming particles: \(\mathrm{fh}\equiv \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 \(\langle \sigma v\rangle _{\mathrm{p}}\) (Eq. (21)) instead of \(\langle \sigma v\rangle _{\mathrm{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. 3 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 \(\phi \) and a Dirac fermion \(\psi \). Both new fields are charged only under a global \(U(1)_{\mathrm{dark}}\) symmetry with charges: \(Q_{\mathrm{dark}}(\psi )=1\), \(Q_{\mathrm{dark}}(\phi )=2\). A non-zero VEV of \(\phi \) spontaneously breaks this symmetry to the discrete \({\mathbb {Z}}_2\) parity. The Goldstone boson associated with this symmetry breaking contributes to DR while the lighter dark fermion (\(\psi \) 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 $$\begin{aligned} \begin{aligned} {\mathcal {L}}_{H,\phi }&= \left( D_\mu H\right) ^\dagger \left( D^\mu H\right) + \mu _H^2H^\dagger H- \lambda _H(H^\dagger H)^2\\&\quad + \partial _\mu \phi ^*\partial ^\mu \phi + \mu _\phi ^2(\phi ^*\phi )^2\\&\quad -\lambda _\phi (\phi ^*\phi )^2 -\kappa (H^\dagger H)(\phi ^*\phi ), \end{aligned} \end{aligned}$$ (22) where H is the Standard Model Higgs doublet. After the spontaneous symmetry breaking the scalar fields may be written as $$\begin{aligned} H=\begin{pmatrix} G^+\\ v_H + \frac{{\tilde{h}}+iG^0}{\sqrt{2}} \end{pmatrix},\quad \quad \phi = v_\phi + \frac{{\tilde{\rho }}+i\sigma }{\sqrt{2}}, \end{aligned}$$ (23) where \(\sigma \) is the Goldstone boson of the broken \(U(1)_{\mathrm{dark}}\) and contributes to DR. Two neutral scalars \({\tilde{h}}\) and \({\tilde{\rho }}\) mix, forming two mass eigenstates (h and \(\rho \)): $$\begin{aligned} h={\tilde{h}}\cos \theta +{\tilde{\rho }}\sin \theta ,\quad \quad \rho = {\tilde{\rho }}\cos \theta -{\tilde{h}}\sin \theta , \end{aligned}$$ (24) with masses $$\begin{aligned} m_h^2= & {} 4\lambda _H v_H^2\cos ^2\theta + 4\lambda _\phi v_\phi ^2\sin ^2\theta + 2\kappa \, v_H v_\phi \sin 2\theta , \nonumber \\\end{aligned}$$ (25) $$\begin{aligned} m_\rho ^2= & {} 4\lambda _\rho v_\rho ^2\cos ^2\theta + 4\lambda _H v_H^2\sin ^2\theta - 2\kappa \, v_H v_\phi \sin 2\theta ,\nonumber \\ \end{aligned}$$ (26) and the mixing angle given by the condition $$\begin{aligned} \tan 2\theta =\frac{\kappa \, v_Hv_\phi }{\lambda _H v_H^2 - \lambda _\phi v_\phi ^2}. \end{aligned}$$ (27) 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, \(\lambda _\phi \) and \(\kappa \), and the mass of the non-SM-like scalar, \(m_\rho \). Other parameters present in the mass formulas (25) and (26) may be expressed as $$\begin{aligned} \sin \theta\simeq & {} \frac{\kappa \, m_\rho v_H}{\left[ \lambda _\phi (m_h^2-m_\rho ^2)^2 - \kappa ^2v_H^2(m_h^2-m_\rho ^2)\right] ^{1/2}}\nonumber \\\approx & {} \frac{\kappa }{\lambda _\phi ^{1/2}}\left( \frac{m_\rho }{90\;\,\mathrm GeV}\right) , \end{aligned}$$ (28) $$\begin{aligned} \lambda _H= & {} \frac{m_h^2\cos ^2\theta + m_\rho ^2\sin ^2\theta }{4v_H^2}\approx 0.13\cos ^2\theta , \end{aligned}$$ (29) $$\begin{aligned} v_\phi ^2= & {} \frac{m_\rho ^2\cos ^2\theta + m_h^2\sin ^2\theta }{4\lambda _\phi }\approx m_\rho ^2\left( \frac{\cos \theta }{2\lambda _\phi ^{1/2}}\right) ^2, \end{aligned}$$ (30) where the approximate equalities are valid if the scalar \(\rho \) is relatively light and does not mix strongly with the SM-like Higgs i.e. \(m^2_\rho \ll m^2_h\), \(\sin ^2\theta \ll 1\). Let us now estimate phenomenologically interesting range of values for \(\kappa \) and \(\lambda _\phi \) that will be used in our numerical scan. The contribution from the dark particles to the invisible Higgs decay width may be approximated as $$\begin{aligned} \Delta \Gamma _{h,\mathrm{inv}}\ge & {} \Gamma (h\rightarrow \sigma \sigma )\nonumber \\&+ \Gamma (h\rightarrow \rho \rho ) \approx \frac{\sin ^2\theta }{32\pi }\frac{m_h^3}{v_\phi ^2} \approx \frac{\kappa ^2}{8\pi \cos ^2\theta }\frac{v_H^2}{m_h},\nonumber \\ \end{aligned}$$ (31) where we conservatively neglected the Higgs decays into dark fermions. The LHC constraint on the invisible Higgs decays i.e. $$\begin{aligned} \Delta \Gamma _{h,\mathrm{inv}}< \frac{B_{\mathrm{inv}}\cos ^2\theta }{1 - B_{\mathrm{inv}}}\Gamma _{h,\mathrm SM}, \end{aligned}$$ (32) where \(B_{\mathrm{inv}}\lesssim 24\%\) [19] and \(\Gamma _{h,\mathrm SM}\simeq 4\) MeV, may be used together with Eq. (31) to obtain the following upper bound on the portal coupling: $$\begin{aligned} \kappa \lesssim 0.01. \end{aligned}$$ (33) Bounds on the singlet scalar self-coupling, \(\lambda _\phi \), are much weaker. They follow from the requirement of perturbativity to some high energy scale and may be obtained from the RG equations [20]: $$\begin{aligned} 16\pi ^2\frac{\mathrm{d}\lambda _\phi }{\mathrm{d}t}&= 10\lambda _\phi ^2 + \kappa ^2,\nonumber \\ 16\pi ^2\frac{\mathrm{d}\kappa }{\mathrm{d}t}&= \kappa \left( 2\kappa + 4\lambda _\phi + 6\lambda _H + 3 y_t^2 -\frac{9}{4}g_2^2 - \frac{3}{4}g_1^2\right) ,\nonumber \\ \end{aligned}$$ (34) where \(t\equiv Q^2/Q_0^2\). The scalar couplings stay perturbative up to the GUT energy scale if at low energy \(\lambda _\phi \lesssim 0.25\). On the other hand, if \(\lambda _\phi =1\) at the weak scale, it becomes non-perturbative already at scale of order 100 TeV. Both values i.e. \(\lambda _\phi = 0.25,\,1\) will be our reference points in numerical scans in the next section. 4 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 \(\sigma \) one needs to know its cross section for the annihilation into SM particles. In the original paper [6] it was assumed that \(\sigma \) decouples at temperature (just) above the muon mass. Then, in some analyses (e.g. [18, 21]) only the annihilation of \(\sigma \) 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 \(\rho \) exchanged in the s channel may be approximated as3 $$\begin{aligned}&\sum _{\mathrm{spins}}|M_{\sigma \sigma \rightarrow \mu ^+\mu ^-}(s)|^2 \approx 2\kappa ^2s^2\nonumber \\&\quad \times \frac{m_{\mu ^\pm }^2(s-4m_{\mu ^\pm }^2)}{\left[ (s-m_\rho ^2)^2 + \Gamma _\rho ^2m_\rho ^2\right] \left[ (s-m_h^2)^2 + \Gamma _h^2m_h^2\right] }, \end{aligned}$$ (35) $$\begin{aligned}&\sum _{\mathrm{spins}}|M_{\sigma \sigma \rightarrow \pi \pi }(s)|^2 \approx 2\kappa ^2s^2\nonumber \\&\quad \times \frac{\frac{1}{27}(s+\frac{11}{2}m_\pi ^2)^2}{\left[ (s-m_\rho ^2)^2 + \Gamma _\rho ^2m_\rho ^2\right] \left[ (s-m_h^2)^2 + \Gamma _h^2m_h^2\right] }. \end{aligned}$$ (36) We checked that the effects from \(\tau ^\pm \), gluons and free quarks become important at temperatures \(T\gtrsim 0.4\div 0.5\) GeV. This suggests that the Goldstone bosons stay in equilibrium at temperatures above \(\Lambda _{\mathrm{QCD}}\). Thus, we assume in our analysis that \(\sigma \) 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\sim m_\rho ^2\). If the resonance is narrow, i.e. when \(\left( {\Gamma _\rho }/{m_\rho }\right) ^2\ll 1\), one may use the approximation $$\begin{aligned} \frac{1}{\pi }\frac{m_\rho \Gamma _\rho }{(s-m_\rho ^2)^2 + \Gamma _\rho ^2m_\rho ^2} \longrightarrow \delta (s-m_\rho ^2). \end{aligned}$$ (37) We assume for simplicity that the CDM fermion is not very light and fulfills the condition \(m_{\mathrm{CDM}}>m_\rho /2\). Then the width of \(\rho \) is dominated by its decays into Goldstone bosons pairs, giving $$\begin{aligned} \Gamma _\rho \approx \Gamma (\rho \rightarrow \sigma \sigma )= \frac{m_\rho ^3\cos ^2\theta }{64\pi v_\phi ^2}\approx m_\rho \left( \frac{\lambda _\phi }{16\pi }\right) . \end{aligned}$$ (38) Because \(\left( {\Gamma _\rho }/{m_\rho }\right) ^2=(\lambda _\phi /16\pi ^2)^2\ll 1\) approximation (37) may be used even for \(\lambda _\phi \) as big as 1. We do not consider larger values of \(\lambda _\phi \) 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 $$\begin{aligned} \sum _{\mathrm{spins}}|M_{\sigma \sigma \rightarrow \mu ^+\mu ^-}|^2= & {} 2\pi \kappa ^2m_\rho ^3\,\delta (s-m_\rho ^2)\nonumber \\&\times \frac{m_{\mu ^\pm }^2(m_\rho ^2-4m_{\mu ^\pm }^2)}{\Gamma _\rho \left[ (m_\rho ^2-m_h^2)^2 + \Gamma _h^2m_h^2\right] }, \end{aligned}$$ (39) $$\begin{aligned} \sum _{\mathrm{spins}}|M_{\sigma \sigma \rightarrow \pi \pi }|^2= & {} 2\pi \kappa ^2m_\rho ^3\,\delta (s-m_\rho ^2)\nonumber \\&\times \frac{\frac{1}{27}(m_\rho ^2+\frac{11}{2}m_\pi ^2)^2}{\Gamma _\rho \left[ (m_\rho ^2-m_h^2)^2 + \Gamma _h^2m_h^2\right] }. \end{aligned}$$ (40) Using this approximation we can also rewrite integrals (11), (20) and (21) as $$\begin{aligned} {\mathscr {D}}\Phi= & {} \frac{1}{2x^2} \int _{1}^\infty \mathrm{d}y\frac{1}{\sqrt{y^2-1}}\; \frac{1}{(|s_{\mathrm{out}}|-e^{-p})(e^{p+2z}-|s_{\mathrm{out}}|)}\nonumber \\&\times \ln \left[ \frac{\cosh \left( \frac{1}{2}(p+q)+z\right) + s_{\mathrm{in}}}{\cosh \left( \frac{1}{2}(p-q)+z\right) + s_{\mathrm{in}}} \right] \nonumber \\&\times \ln \left[ \frac{\cosh \left( \frac{1}{2}\big (p+qV\big )\right) + s_{\mathrm{out}}}{\cosh \left( \frac{1}{2}\big (p-qV\big )\right) + s_{\mathrm{out}}} \right] , \end{aligned}$$ (41) $$\begin{aligned} \langle \sigma v\rangle _{\mathrm{MB}}= & {} \frac{\sum _{\mathrm{spins}}|M|^2}{512\pi }\frac{x^5}{m_{\mu ^\pm }^2} \sqrt{\frac{m_\rho ^2}{m_{\mu ^\pm }^2}-4} \;\mathrm{K}_1\left( \frac{x\,m_\rho }{m_{\mu ^\pm }}\right) , \end{aligned}$$ (42) $$\begin{aligned} \langle \sigma v\rangle _{\mathrm{p}}= & {} \frac{\sum _{\mathrm{spins}}|M|^2}{512\pi }\frac{x^5}{m_{\mu ^\pm }^3}\frac{1}{(\zeta ^{\mathrm{in}})^2} \sqrt{1-\frac{4m_{\mu ^\pm }^2}{m_\rho ^2}} \nonumber \\&\times \int _{m_\rho }^\infty \mathrm{d}E_+\frac{e^{-\frac{x}{2m_{\mu ^\pm }}E_+}}{\sinh \left( \frac{x}{2m_{\mu ^\pm }}E_+\right) }\nonumber \\&\times \ln \left[ \frac{\mathrm{fh}\left( \frac{x}{4m_{\mu ^\pm }} \left( E_++\sqrt{E_+^2-m_\rho ^2}\right) \right) }{\mathrm{fh}\left( \frac{x}{4m_{\mu ^\pm }}\left( E_+-\sqrt{E_+^2-m_\rho ^2}\right) \right) } \right] , \end{aligned}$$ (43) 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_{\mathrm{eff}}\). It is usually defined by the formula $$\begin{aligned} \rho _R = \left[ 1 + \left( \frac{\rho _\nu ^{\mathrm{ifo}}}{\rho _\gamma ^{\mathrm{ifo}}}\right) N_{\mathrm{eff}}\right] \rho _\gamma = \left[ 1 + \frac{7}{8}\left( \frac{4}{11}\right) ^{4/3} N_{\mathrm{eff}}\right] \rho _\gamma ,\nonumber \\ \end{aligned}$$ (44) where \(\rho _R\), \(\rho _\gamma \) and \(\rho _\nu \) correspond to the present total, photons and neutrinos radiation energy density, respectively. Index ifo denotes instantaneous freeze-out approximation. In general we can write $$\begin{aligned} \rho _R = \rho _\gamma + \sum _{i=1}^3\rho _{\nu _i} + \rho _X, \end{aligned}$$ (45) where \(\rho _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=\sigma \). \(N_{\mathrm{eff}}\) is normalized in such a way that using the instantaneous freeze-out approximation for neutrino decoupling and putting \(\rho _X=0\) we get \(N_{\mathrm{eff}}=3\). In the Standard Model due to the lack of perfect kinetic equilibrium during neutrino decoupling one gets \(N_{\mathrm{eff}}^{\mathrm{SM}}\approx 3.045\) [25]. For simplicity, we will define additional effective number of neutrino species as $$\begin{aligned} \Delta N_{\mathrm{eff}}\equiv N_{\mathrm{eff}} - N_{\mathrm{eff}}^{\mathrm{SM}}, \end{aligned}$$ (46) which for the Goldstone boson \(\sigma \) (additional factor \(\frac{7}{8}\) due to different statistics and number of degrees of freedom as compared to \(\nu \)) gives $$\begin{aligned} \Delta N_{\mathrm{eff}}= & {} \frac{4}{7}\left( \frac{T_\gamma ^{\mathrm{ifo}}}{T_\nu ^{\mathrm{ifo}}}\right) ^4\frac{\rho _\sigma }{\rho _\gamma } = \left. \frac{4}{7}\frac{\rho _\sigma }{\rho _\gamma }\right| _{T=T_{\mathrm{end}}} \nonumber \\= & {} \left. \frac{4}{7}\left( \frac{Y}{Y_{\mathrm{eq}}}\right) ^{\frac{4}{3}}\right| _{x=x_{\mathrm{end}}}, \end{aligned}$$ (47) where \(x_{\mathrm{end}}\approx 20\) corresponds to \(T_{\mathrm{end}}\approx 5\) MeV (temperature before neutrino decoupling when the freeze-out process of \(\sigma \) is completed). 4.1 Instantaneous freeze-out approximation Instantaneous freeze-out approximation is one of the simplest methods for relic density estimation. By definition, freeze-out takes place at such \(x_f\) for which the freeze-out parameter \(\eta (x)\), being the ratio between the interaction rate \(\Gamma \) and the Hubble parameter H, equals one: $$\begin{aligned} \eta (x) = \left. \frac{\Gamma }{H}\right| _x = \left. \frac{n\langle \sigma v\rangle }{H}\right| _x,\quad \quad \eta (x_f)=1. \end{aligned}$$ (48) At temperatures below \(T_f\) (\(x>x_f\)) the yield of considered particles remains constant i.e. \(Y(x)=Y^{\mathrm{eq}}(x_f)\). For \(\langle \sigma v\rangle = \langle \sigma v\rangle _{\mathrm{MB}}\) (see Eq. (20)) with matrix elements squared as in (39)–(40) and \(\Gamma _\rho \) given by (38), one gets the following condition for \(x_f\) Open image in new window Fig. 1 Top: contour lines of constant \(\Delta N_{\mathrm{eff}}\) in (\(m_\rho \), \(\kappa \)) plane for \(\lambda _\phi =0.25\) (left panel) and \(\lambda _\phi =1\) (right panel) obtained with Eq. (48) in Maxwell-Boltzmann approximation (dark/light blue dashed lines include annihilation into muons and pions/muons only) and Eq. (12) with full particles statistics (red lines – see details in Sect. 4.2). Black lines correspond to astrophysical bounds [26]: for Goldstone boson free-streaming out of the proto-neutron stars (A) and for energy loss due to exotic species in SN 1987A in two different approximations (B). Green lines denote current (LHC – Eq. (33)) and future (ILC for \(B_{\mathrm{inv}}=0.3\%\)) limit on invisible Higgs decay. We did not include effects from the CDM sector. Bottom: \(\eta (T)\) dependence (see Eq. (49)) for large dots marked in the upper panels – brown (black) lines correspond to \(\kappa =10^{-3}\; (2\cdot 10^{-4})\), whereas solid (dashed) to \(m_\rho =0.4\; (1)\) GeV. The freeze-out process for parameters corresponding to the black dots for \(\kappa =2\cdot 10^{-4}\) is presented in more detail in Fig. 2 $$\begin{aligned} {\mathcal {C}} \frac{x_f^4}{\sqrt{g_*(x_f)}} \mathrm{K}_1\left( \frac{x_f\, m_\rho }{m_{\mu ^\pm }}\right) =1, \end{aligned}$$ (49) where $$\begin{aligned} {\mathcal {C}}\equiv & {} \frac{\kappa ^2}{\lambda _\phi } \frac{\sqrt{45}\,\zeta (3)}{32 \pi ^{5/2}} \frac{m_\rho ^5 M_{\mathrm{Pl}}}{m_h^4 m_{\mu ^\pm }^2} \nonumber \\&\times \left[ \left( 1-\frac{4m_{\mu ^\pm }^2}{m_\rho ^2}\right) ^{3/2}+ \frac{1}{27} \frac{m_\rho ^2}{m_{\mu ^\pm }^2} \left( 1+\frac{11}{2}\frac{m_\pi ^2}{m_\rho ^2}\right) ^{2} \left( 1-\frac{4m_{\mu ^\pm }^2}{m_\rho ^2}\right) ^{1/2} \right] .\nonumber \\ \end{aligned}$$ (50) 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 \(\Delta N_{\mathrm{eff}}\) obtained in the instantaneous freeze-out approximation, \(Y=Y_{\mathrm{eq}}(x_f)\), as a function of \(m_\rho \) and \(\kappa \) (blue lines), using definition (47) and the freeze-out temperature calculated from Eq. (49). The results are compared to those obtained by solving the Boltzmann equation (red lines) – see Sect. 4.2. In almost all cases the instantaneous freeze-out approximation overestimates \(\Delta N_{\mathrm{eff}}\) for a given values of \(\kappa \) and \(\lambda _\phi \) (dark blue lines are below corresponding red lines). The best agreement between the instantaneous freeze-out approximation and the full result is achieved for \(\Delta N_{\mathrm{eff}}\) close to 0.1, although in that region the instantaneous freeze-out approximation works only for \(m_\rho \gtrsim 1\) GeV. Note also that neglecting the \(\sigma \) 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 \(\Delta N_{\mathrm{eff}}\) (the difference amplifies when \(\kappa \) decreases).5 Open image in new window Fig. 2 Yield Y (upper panels) and pseudopotential z (lower panels) of Goldstone boson \(\sigma \) as functions of temperature for \(\lambda _\phi =0.25\), \(\kappa =2\cdot 10^{-4}\) and for two values of \(m_\rho \) – see black dots in the upper left panel of Fig. 1. Line colors (corresponding to the methods used) are described in the legends. Dashed blue line in the upper right panel denotes the instantaneous freeze-out approximation result (see Fig. 1) It is worth emphasizing that the instantaneous freeze-out approximation breaks down for very small \(\kappa \) and \(m_\rho \lesssim 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 \(\kappa \) the \(\eta (T)\) function never reaches 1. The current experimental limits, shown in Fig. 1, exclude \(\kappa \gtrsim 10^{-2}\) and the region of small \(m_\rho \) with \(\kappa \gtrsim 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]. 4.2 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 \(\kappa \), where the last approximation breaks down, are especially important from the viewpoint of near future experiments that will be able to probe \(\Delta N_{\mathrm{eff}}\) with accuracy even better than 0.03. Let us start the discussion by considering two sample points from the parameter space with small \(\kappa \): 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\sim 100\div 200\) MeV and \(30\div 50\) MeV is characterized by strong variability. As a result, it is harder for \(\sigma \) particles to follow the equilibrium density i.e. larger \(\kappa \) 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 \(\sigma \) freeze-out changes relatively mildly (e.g. for \(T\gtrsim 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_\sigma (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 \(\sigma \) freezes in for \(T\sim 60\div 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_\rho \) reaches its maximum in this temperature range (see also \(\eta (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 as6 $$\begin{aligned} {\mathscr {D}}\Phi _{\mathrm{fBE}}\approx & {} {\mathscr {D}}\Phi _{\mathrm{MB}}\times \zeta ^{\mathrm{in}}, \nonumber \\ {\mathscr {D}}\Phi _{\mathrm{fFD}}\approx & {} {\mathscr {D}}\Phi _{\mathrm{MB}}\times \zeta ^{\mathrm{in}}, \end{aligned}$$ (51) $$\begin{aligned} {\mathscr {D}}\Phi _{\mathrm{pBE}}\approx & {} {\mathscr {D}}\Phi _{\mathrm{MB}}\times \frac{1}{\zeta ^{\mathrm{in}}} \coth \left( \frac{m_\rho x}{4 m_{\mu ^\pm }}\right) \nonumber \\ {\mathscr {D}}\Phi _{\mathrm{pFD}}\approx & {} {\mathscr {D}}\Phi _{\mathrm{MB}}\times \frac{1}{\zeta ^{\mathrm{in}}} \tanh \left( \frac{m_\rho x}{4 m_{\mu ^\pm }}\right) , \end{aligned}$$ (52) $$\begin{aligned} {\mathscr {D}}\Phi _{\mathrm{BEFD}}\approx & {} {\mathscr {D}}\Phi _{\mathrm{MB}}\nonumber \\ {\mathscr {D}}\Phi _{\mathrm{FDBE}}\approx & {} {\mathscr {D}}\Phi _{\mathrm{MB}}, \end{aligned}$$ (53) $$\begin{aligned} {\mathscr {D}}\Phi _{\mathrm{BEBE}}\approx & {} {\mathscr {D}}\Phi _{\mathrm{MB}}\times \coth ^2\left( \frac{m_\rho x}{4 m_{\mu ^\pm }}\right) \nonumber \\ {\mathscr {D}}\Phi _{\mathrm{FDFD}}\approx & {} {\mathscr {D}}\Phi _{\mathrm{MB}}\times \tanh ^2\left( \frac{m_\rho x}{4 m_{\mu ^\pm }}\right) , \end{aligned}$$ (54) where (see Eq. (42)) $$\begin{aligned} {\mathscr {D}}\Phi _{\mathrm{MB}}= e^{-2z}\frac{\sum _{\mathrm{spins}}|M|^2}{2x} \sqrt{\frac{m_\rho ^2}{m_{\mu ^\pm }^2}-4}\; \mathrm{K}_1\left( \frac{x\,m_\rho }{m_{\mu ^\pm }}\right) . \end{aligned}$$ (55) 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\left( m_\rho x/4 m_{\mu ^\pm }\right) \)) starts to be important for small \(m_\rho \) and x (larger T) and also for small \(\kappa \). The dependence on \(\kappa \) follows from the fact that \(\sigma \) with smaller \(\kappa \) decouples at higher temperature i.e. at smaller x. These effects can be seen in Fig. 3, where for small enough \(m_\rho \) and \(\kappa \) red curves are placed partially under the blue ones i.e. in order to get a given value of \(\Delta N_{\mathrm{eff}}\) smaller \(\kappa \) is needed when using the Boltzmann equation with full statistics effects as compared to the MB approximation. Open image in new window Fig. 3 As top panels in Fig. 1 but for solutions of the Boltzmann equation with and without approximations defined in Sect. 2.1. Colors for continuous lines are described in Fig. 2. Red dashed lines correspond to calculation with the pions contribution neglected One can also observe from Fig. 3 that for a given value of \(\kappa \) and \(m_\rho \) in most cases \(\Delta N_{\mathrm{eff}}^{\mathrm{pBE}}<\Delta N_{\mathrm{eff}}^{\mathrm{MB}}<\Delta N_{\mathrm{eff}}^{\mathrm{fBE}}\) and relative differences between approximations decrease when \(\kappa \) increases. Only for small enough \(\kappa \) and \(m_\rho \) we have \(\Delta N_{\mathrm{eff}}^{\mathrm{pBE}}>\Delta N_{\mathrm{eff}}^{\mathrm{MB}}\). 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_\rho \) and \(\kappa \), where pions statistics matters) full inclusion of incoming and outgoing statistics results in smaller \(\Delta N_{\mathrm{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_{\mathrm{eq}}(x)\)) which causes weaker change in \(|\mathrm{d}z/\mathrm{d}x|\) 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 by \({\tilde{S}}_I(z,x)\)) leading to faster \(|\mathrm{d}z/\mathrm{d}x|\) 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_{\mathrm{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 \(\sigma \) 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_\rho \le 2m_\pi \) (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 \(\Delta N_{\mathrm{eff}}\sim 0.05\) and more.7 Thus, statistics of both incoming and outgoing particles are relevant, especially for moderate values of \(\kappa \) and \(m_\rho \). Minimal value of \(\Delta N_{\mathrm{eff}}\) obtained in the scan (\(\sim 0.06\)) is related to the starting moment for the calculation i.e. \(T_{\mathrm{start}}=400\) MeV, however for larger T \(Y_{\mathrm{eq}}(T)\) does not change significantly. Therefore, if \(\sigma \) freezes out during or after the QCD phase transition, it shall give contribution to \(N_{\mathrm{eff}}\) measurable by near future experiments. On the other hand, maximal \(\Delta N_{\mathrm{eff}}\) equals approximately 0.5 and is achieved for large \(\kappa \) and \(m_\rho \) 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 \(\Delta N_{\mathrm{eff}}\sim 0.39\) for moderate values of \(\lambda _\phi \) 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 \(\lambda _\phi \) may be more challenging in this context as the lines of constant \(\Delta N_{\mathrm{eff}}\) move towards smaller \(\kappa \), 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 for8 \(m_\rho \le 1.2\) GeV, \(\kappa \le 10^{-3}\) and taking into account only (dominant) annihilation into pions (Eq. (40)). Open image in new window Fig. 4 Contour lines of constant \(\Delta N_{\mathrm{eff}}\) (top lines: 0.3, bottom lines: 0.1) in (\(m_\rho \), \(\kappa \)) plane for \(\lambda _\phi =0.25\) and bosons/fermions annihilation (left/right panel). Continuous (dashed) red lines correspond to mixed (homogeneous) statistics – see Eqs. (53) and (54) e.g. BEBE case refers to \(\sigma \) annihilation into pions in Weinberg’s Higgs portal model. Colors for other lines (MB, fBE/fFD, pBE/pFD) are the same as in the previous plots 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 \({\mathscr {D}}\Phi \) (see Eq. (53)). The main reason for smaller/bigger values of \(\Delta N_{\mathrm{eff}}\) in such mixed cases, especially when compared to the results of the MB approximation, is the presence of the statistics factor \(s_{\mathrm{in}}={\pm }\,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\left( m_\rho x/4 m_{\mu ^\pm }\right) \) and \(\tanh ^2\left( m_\rho x/4 m_{\mu ^\pm }\right) \)) in \({\mathscr {D}}\Phi \) given by (54)) are important for small exchanged particle mass (here \(m_\rho \)) and coupling (here \(\kappa \)). The effects of outgoing particle statistics decrease with \(m_\rho \) and \(\kappa \) and eventually become negligible. 5 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 \(\Delta N_{\mathrm{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_\rho \) and \(\kappa \) for which the Boltzmann equation must be used. In most of the remaining regions of the parameter space the instantaneous freeze-out approximation overestimates \(\Delta N_{\mathrm{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 \(\Delta N_{\mathrm{eff}}\sim 0.1\) and \(m_\rho \gtrsim 1\) GeV. The resonant exchange of the scalar \(\rho \) (for \(m_\rho \) up to a few hundreds MeV) also plays an important role. When calculating the annihilation cross section of DR particles (Goldstone bosons \(\sigma \)) it is crucial to include muons and pions as the final states. Pions are quite often ignored which may lead to underestimation of \(\Delta N_{\mathrm{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 \(\Delta N_{\mathrm{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 \(\Delta N_{\mathrm{eff}}\) than inclusion of only some of the effects – those in evaluation of \(Y_{\mathrm{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 \(\Delta N_{\mathrm{eff}}\). One may expect that results of near future experiments will lead to much better determination of \(\Delta N_{\mathrm{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 \(\rho \) when DR particles \(\sigma \) may decouple kinetically before thermal decoupling. This might happen because in elastic scattering particle \(\rho \) 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. Footnotes 1. Note the difference with respect to \(J_n\) function defined in Eq. (10) in [15]. 2. If \(\chi \) and \({\bar{\chi }}\) are distinguishable one should consider \(n = n_\chi + n_{{\bar{\chi }}}\) instead. 3. We neglected the mass difference between the charged and neutral pions and summed effectively their contributions at the level of amplitudes instead of cross sections. 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_\rho /m_{\mu ^\pm })y\), \(q=x(m_\rho /m_{\mu ^\pm })\sqrt{y^2-1}\) and \(V=\sqrt{1-4(m_{\mu ^\pm }/m_\rho )^2}\). 5. One can see that for larger \(\Delta N_{\mathrm{eff}}\) the freeze-out approximation gives better agreement with the Boltzmann equation solutions when only \(\sigma \) annihilation into \(\mu ^\pm \) is taken into account, than including also annihilation into \(\pi ^\pm \) and \(\pi ^0\). In particular, for \(\Delta N_{\mathrm{eff}}=0.3\) one can observe a quite accurate coincidence in almost whole range of \(m_\rho \). But this is just an example of situations when effects resulting from two wrong assumptions/approximations partially cancel out, leading to a final result being close to the correct one. 6. In four letter subscripts (BEFD, FDBE, BEBE and FDFD) the two first (last) letters denote the statistics of incoming (outgoing) particles. 7. This effect for a Goldstone boson with just one degree of freedom is of similar magnitude as the summed effect of kinetic non equilibrium during decoupling for six neutrino degrees of freedom in the SM (\(N_{\mathrm{eff}}^{\mathrm{SM}}\approx 3.045\)). 8. Asymptotic behavior for larger \(m_\rho \) and \(\kappa \) resembles that from the previous plots. Notes Acknowledgements This work has been partially supported by National Science Centre, Poland, under research Grants DEC-2014/15/B/ST2/02157 and DEC-2012/04/A/ST2/00099. MO acknowledges partial support from National Science Centre, Poland, grant DEC-2016/23/G/ST2/04301. PS acknowledges support from National Science Centre, Poland, grant DEC-2015/19/N/ST2/01697. References 1. P.A.R. Ade et al., [Planck Collaboration]. Astron. Astrophys. 594, A13 (2016). arXiv:1502.01589 [astro-ph.CO] 2. A. G. Riess et al., Astrophys. J. 826(1), 56 (2016). arXiv:1604.01424 [astro-ph.CO] 3. R.C. Nunes, A. Bonilla, Mon. Not. Roy. Astron. Soc. 473, 4404 (2018). arXiv:1710.10264 [astro-ph.CO]ADSCrossRefGoogle Scholar 4. A. Heavens, Y. Fantaye, E. Sellentin, H. Eggers, Z. Hosenie, S. Kroon, A. Mootoovaloo, Phys. Rev. Lett. 119(10), 101301 (2017). arXiv:1704.03467 ADSCrossRefGoogle Scholar 5. K. N. Abazajian et al. [CMB-S4 Collaboration]. arXiv:1610.02743 [astro-ph.CO] 6. S. Weinberg, Phys. Rev. Lett 110(24), 241301 (2013). arXiv:1305.1971 [astro-ph.CO]ADSCrossRefGoogle Scholar 7. S. Hannestad, J. Madsen, Phys. Rev. D 52, 1764 (1995). arXiv:astro-ph/9506015 ADSCrossRefGoogle Scholar 8. A.D. Dolgov, S.H. Hansen, D.V. Semikoz, Nucl. Phys. B 503, 426 (1997). arXiv:hep-ph/9703315 ADSCrossRefGoogle Scholar 9. N.Y. Gnedin, O.Y. Gnedin, Astrophys. J. 509, 11 (1998). arXiv:astro-ph/9712199 ADSCrossRefGoogle Scholar 10. G. Mangano, G. Miele, S. Pastor, T. Pinto, O. Pisanti, P.D. Serpico, Nucl. Phys. B 729, 221 (2005). arXiv:hep-ph/0506164 ADSCrossRefGoogle Scholar 11. S. Esposito, G. Miele, S. Pastor, M. Peloso, O. Pisanti, Nucl. Phys. B 590, 539 (2000). arXiv:astro-ph/0005573 ADSCrossRefGoogle Scholar 12. G. Mangano, G. Miele, S. Pastor, M. Peloso, Phys. Lett. B 534, 8 (2002). arXiv:astro-ph/0111408 ADSCrossRefGoogle Scholar 13. J. Birrell, J. Wilkening, J. Rafelski, J. Comput. Phys 281, 896 (2014). arXiv:1403.2019 [math.NA]ADSCrossRefGoogle Scholar 14. J. Birrell, C.T. Yang, P. Chen, J. Rafelski, Phys. Rev. D 89, 023008 (2014). arXiv:1212.6943 [astro-ph.CO]ADSCrossRefGoogle Scholar 15. A.D. Dolgov, K. Kainulainen, Nucl. Phys. B 402, 349 (1993). arXiv:hep-ph/9211231 ADSCrossRefGoogle Scholar 16. J. Bernstein, L.S. Brown, G. Feinberg, Phys. Rev. D 32, 3261 (1985)ADSCrossRefGoogle Scholar 17. P. Gondolo, G. Gelmini, Nucl. Phys. B 360, 145 (1991)ADSCrossRefGoogle Scholar 18. C. Garcia-Cely, A. Ibarra, E. Molinaro, JCAP 1311, 061 (2013). arXiv:1310.6256 [hep-ph]ADSCrossRefGoogle Scholar 19. V. Khachatryan et al., [CMS Collaboration]. JHEP 1702, 135 (2017). arXiv:1610.09218 [hep-ex] 20. W .F. Chang, J .N. Ng, JCAP 1607(07), 027 (2016). arXiv:1604.02017 [hep-ph]ADSCrossRefGoogle Scholar 21. L .A. Anchordoqui, P .B. Denton, H. Goldberg, T .C. Paul, L H M Da Silva, B .J. Vlcek, T .J. Weiler, Phys. Rev. D 89(8), 083513 (2014). arXiv:1312.2547 [hep-ph]ADSCrossRefGoogle Scholar 22. M. B. Voloshin, Sov. J. Nucl. Phys. 44 (1986) 478 [Yad. Fiz. 44 (1986) 738]Google Scholar 23. J.F. Gunion, H.E. Haber, G.L. Kane, S. Dawson, Front. Phys. 80, 1 (2000)Google Scholar 24. K .W. Ng, H. Tu, T .C. Yuan, JCAP 1409(09), 035 (2014). arXiv:1406.1993 [hep-ph]ADSCrossRefGoogle Scholar 25. P. F. de Salas and S. Pastor, JCAP 1607(07), 051 (2016). arXiv:1606.06986 [hep-ph] 26. H. Tu, K.W. Ng, JHEP 1707, 108 (2017). arXiv:1706.08340 [hep-ph]ADSCrossRefGoogle Scholar 27. M. Drees, F. Hajkarim and E. R. Schmitz, JCAP 1506(06), 025 (2015). arXiv:1503.03513 [hep-ph] 28. T. Binder, T. Bringmann, M. Gustafsson and A. Hryczuk, Phys. Rev. D 96(11), 115010 (2017). arXiv:1706.07433 [astro-ph.CO] Copyright information © The Author(s) 2018 Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP3 Authors and Affiliations Marek Olechowski1Paweł  Szczerbiak1Email author1.Institute of Theoretical Physics, Faculty of PhysicsUniversity of WarsawWarsawPoland


This is a preview of a remote PDF: https://link.springer.com/content/pdf/10.1140%2Fepjc%2Fs10052-018-6195-0.pdf

Marek Olechowski, Paweł Szczerbiak. Effects of quantum statistics on relic density of dark radiation, The European Physical Journal C, 2018, 704, DOI: 10.1140/epjc/s10052-018-6195-0