Simple estimate of BBN sensitivity to light freeze-in dark matter

We provide a simple analysis of the big-bang nucleosynthesis (BBN) sensitivity to the light dark matter (DM) generated by the thermal freeze-in mechanism. It is shown that the ratio of the effective neutrino number shift $\Delta N_{\nu}$ over the DM relic density $\omega\equiv \Omega h^2$, denoted by $R_\chi\equiv\Delta N_\nu/\omega$, cancels the decaying particle mass and the feeble coupling, rendering therefore a simple visualization of $\Delta N_{\nu}$ at the BBN epoch in terms of the DM mass. This property drives one to conclude that the shift with a sensitivity of $\Delta N_{\nu}\simeq \mathcal{O}(0.1)$ cannot originate from a single warm DM under the Lyman-$\alpha$ forest constraints. For the cold-plus-warm DM scenarios where the Lyman-$\alpha$ constraints are diluted, the ratio $R_\chi$ can be potentially used to test the thermal freeze-in mechanism in generating a small warm component of DM and a possible excess at the level of $\Delta N_{\nu}\simeq \mathcal{O}(0.01)$.


I. INTRODUCTION
Hidden light particles with typical masses around keV scale are widely considered as the dark matter (DM) candidates [1]. Such light particles are usually associated with feeble couplings and have a non-thermal history in the early Universe. Feeble interactions with the Standard Model (SM) sector are extremely difficult to probe directly at colliders and indirectly in flavor physics. Nevertheless, the precision era of cosmology now provides us with striking opportunities to search indirectly for these hidden particles. Among other things, the precise measurements of light elements produced at the big-bang nucleosynthesis (BBN) epoch set constraints on any new physics effects that can cause a significant shift of the effective neutrino number ∆N ν ≡ N ν − N SM ν away from the SM prediction N SM ν = 3. When the BBN measurements, including the primordial abundances of helium-4 and deuterium, are combined with the cosmic microwave background (CMB) data (CMB+BBN+Y p +D), the constraints become even tighter and N ν is found to be [2]: which sets a 2σ upper bound on the shift ∆N ν < 0.151. The generation of an observable shift of ∆N ν at the BBN epoch or of ∆N ef f at the CMB regime with respect to the SM prediction N SM ef f = 3.045 obtained after taking into account the non-instantaneous neutrino decoupling effects [3][4][5][6][7][8], has opened up an avenue to probe light hidden species in the early Universe (see, e.g., Refs. [9][10][11][12][13][14][15][16]). In particular, for the DM scenarios with a nonthermal history, such as the representative freeze-in evolution [17,18], if the hidden species reside at keV or sub-keV scale, they could contribute to the Hubble expansion as extra radiation prior to the BBN epoch and become non-relativistic before the recombination. Therefore, the BBN measurements can provide a promising avenue to probe and/or constrain the non-thermal light DM which nevertheless has a feeble connection to the visible world.
Here, we present a simple analysis of the BBN sensitivity to the non-thermal light species which is generated by the thermal freeze-in mechanism. While the BBN constraints on non-thermal keV warm DM have been partially discussed by specifying the DM momentum distribution [19,20], we will not refrain ourselves from a single warm DM candidate and extend the analysis to any light hidden species that are difficult to probe directly at terrestrial colliders and indirectly in flavor physics. Noticeably, our analysis could be powerful in testing the mixed cold-plus-warm DM scenarios (see e.g., Refs. [21][22][23][24][25]) in which the small-scale problems encountered in cold DM simulations can be alleviated [26]. As a simple application in this context, our analysis can be used as an economic criterion to test whether a small warm DM component and an observable ∆N ν excess have the common source from a sub-keV hidden species. The analysis presented here is based on the integrated Boltzmann equations of the energy and the particle-number density evolution, instead of solving the unintegrated kinetic equations of momentum distribution.
As will be shown later, the current observation of the DM relic density allows us to see directly how large a shift of ∆N ν can be generated. Remarkably, this visualization is basically independent of the decaying particle mass and the feeble coupling, but essentially dependent on the relic DM mass. Such a freeze-in independence property can promote a simple avenue to test the thermal freeze-in mechanism as the common source of ∆N ν excess and the DM relic density. Given the severe bounds from the Lyman-α forest data (see e.g., Refs. [27][28][29][30][31][32]), we can readily conclude that a single warm DM cannot generate an observable shift of ∆N ν under the current and future sensitivities. Nevertheless, since the Lyman-α constraints are diluted in the mixed cold-pluswarm DM scenarios [22,23,25], we find that the forecast ∆N ν ∼ O(0.01) from a combination of BBN and CMB-S4 [33,34] is able to test a sub-keV DM that makes up only a small fraction of the current relic density.
It should be noted that our analysis presented here could also be generalized to other DM production mechanisms, though it may not be always true to obtain the simple freeze-in independence. Some discussions about the validity of such a freeze-in independence will be presented before the conclusion section.

II. EFFECTIVE NEUTRINO NUMBER SHIFT FROM EXTRA RADIATION
The primary epoch of BBN takes place when the temperature of the early Universe cools down to a few MeV and lasts till O(10) keV (see Ref. [35] for a review). Any extra radiation contributing to the SM plasma at this epoch can be parameterized by the variation of effective neutrino number: which is normalized to the energy density of one active neutrino flavor ρ νi . The main effect of ∆N ν is to modify the Hubble rate through the total energy density, which further affects the freezing temperature of neutron. As a consequence, the neutron-to-proton ratio is altered and a different abundance of helium-4 from what is predicted by the SM with N SM ν = 3 would be generated. In the SM where N ν is fixed, the mass fraction of helium-4, denoted by Y p , depends on the baryon density (η b ) and the neutron mean lifetime (τ n ). When N ν varies due to the extra radiation contribution, one can use the BBN code to obtain a semi-analytic formula of Y p (N ν ) by fixing the two free parameters η b and τ n to their experimental measurements [2,35,36]. Then, a bound on ∆N ν can be extracted by comparing Y p (N ν ) with the measured abundance. With this method, it has been shown that the variation of the neutrino number is currently bounded to be ∆N ν = N ν − 3 < 0.151 at 2σ level from CMB+BBN+Y p +D [2]. However, it should be mentioned that the direct application of this limit to Eq. (2) is valid only for the extra radiation that decouples from the SM plasma before the BBN epoch. For the models, e.g., in which the hidden sector thermalizes with photons, electrons or neutrinos during the BBN process (see, e.g., Refs. [12,13,37]), a simple parameterization of ∆N ν is still absent [38]. In these situations, one should perform a fully numerical simulation via the BBN code (see, e.g., Refs. [35,[39][40][41][42]). For the non-thermal light species considered here, we will assume that they decouple already from the SM plasma before the BBN begins. This allows us to find a simple parametrization of ∆N ν and probe the corresponding signals and/or constraints of the hidden sector. Note that we do not consider the shift of ∆N ef f at the CMB epoch, since we are interested in the light DM that is relativistic at the BBN epoch but becomes non-relativistic prior to the recombination.
In practice, we can express ρ νi in Eq. (2) in terms of the total relativistic energy density of the SM plasma ρ SM : where g ρ * denotes the effective number of relativistic degrees of freedom (d.o.f) for the energy density. Using the fact that the SM plasma conserves entropy after the hidden sector decouples, we can determine the shift ∆N ν in terms of the SM and hidden energy densities: where T ref denotes any reference temperature after the hidden radiation freezes in at T f i and before the species transits to non-relativistic matter at T nr . Here the effective d.o.f for the SM entropy density, g s * (T γeν ) = 10.75, corresponds to the epoch when the relativistic SM plasma contains photons, electrons, positrons and neutrinos. Using Eq. (4), we can then probe and/or constrain new physics effects encoded in ρ rad , by comparing the resulting shift ∆N ν with the bound extracted from BBN precision measurements.

III. PARTICLE-NUMBER AND ENERGY DENSITY EVOLUTION
To see how the freeze-in hidden light species contributes to the shift ∆N ν , we proceed to calculate the particle-number and energy density evolution in light of the integrated Boltzmann equations. For definiteness, we assume that the production is via a thermal scalar decay into two chiral fermions, φ → χ + ψ, from the interaction Lagrangian given by where y 1 is a feeble Yukawa-like coupling. For simplicity, the fermion χ is assumed to be the DM candidate while φ and ψ are non-DM species.
Let us now concentrate on the χ evolution. It is, as usual, convenient to simplify the Boltzmann equation via the particle yield Y χ,n ≡ n χ /s, with the SM entropy density given by The reduced Boltzmann equation for the χ species reads where the Hubble rate is given by H ≈ 1.66 g ρ * T 2 /M Pl , with the Planck mass M Pl ≈ 1.22 × 10 19 GeV. The collision term C χ,n is calculated in the center-of-mass frame of the thermal scalar φ with internal d.o.f g φ = 1, and reads where the Boltzmann distribution f φ = e −E φ /T has been used, and the inverse decay χ + ψ → φ as well as the Pauli-blocking effects have been neglected, with Here K 1 is the first modified Bessel function of the second kind. Integrating Eq. (7) over the temperature T gives the final result of the yield Y χ,n : where we have made the approximation g ρ * ≈ g s * ≈ g * ,χdec , corresponding to the effective d.o.f at the χ decoupling temperature T χdec . The reason behind such an approximation is that, since the yield accumulation is culminated around T χdec O(m φ ), for m φ MeV concerned here, the variations of g ρ * and g s * with respect to the temperature are sufficiently small [43], and the heavier m φ is, the better the approximation will be. In addition, we have fixed the lower limit of the temperature integration to be zero, while a lower limit fixed at the freeze-in temperature T f i m φ results in a negligible numeric difference. This is a generic feature of the freeze-in mechanism as the production is strongly Boltzmann suppressed when T m φ . The evolution of the energy density can be simplified analogously. To this end, we can define the energy yield Y χ,ρ ≡ ρ χ /s 4/3 , so that after χ decouples relativistically from the SM plasma the yield is scale-independent until χ becomes non-relativistic. We have then a Boltzmann equation similar to Eq. (7): where the collision term C χ,ρ in the center-of-mass frame of φ is evaluated to be Note that the approximations made in deriving C χ,n are again applied here. The final result of the energy yield Y χ,ρ is given by It can be immediately seen from Eqs. (9) and (12) that which is now basically independent of the decaying particle mass m φ and the feeble Yukawa-like coupling y (As the decoupling temperature T χdec O(m φ ), an indirectly weak dependence is still observed through the effective d.o.f g * ,χdec , which nevertheless causes only a negligible effect for m φ MeV considered here.). Noticeably, we have checked that this freeze-in independence is also valid when the Bose-Einstein statistics is used for φ (but without the model-dependent thermal mass corrections), up to a small numerical change of the prefactor in Eq. (13). Thus, we can infer from this property that the ratio ∆N ν /ω χ (to be defined later) would essentially depend only on the DM mass m χ . As will be shown later, this renders a fast estimate of how large a shift of ∆N ν can be generated by varying the DM mass within the current relic density ω DM ≡ Ω DM h 2 = 0.120 ± 0.001 [44].

IV. FAST ESTIMATE OF NEUTRINO NUMBER SHIFT
Concerning the shift ∆N ν just prior to the BBN epoch, we can apply Eq. (4) with a reference temperature T ref = T γeν , leading to The relic density of χ at present day is estimated to be Using the current densities s 0 = 2891.2 cm −3 and ρ c = 1.05×10 −5 h 2 ·GeV·cm −3 [45] as input, we can determine the ratio R χ ≡ ∆N ν /ω χ , with the final result given by which depends on the DM mass m χ , and has only a weak dependence on the decaying mass m φ through the effective d.o.f g * ,χdec . Equation (16) is the main result obtained in this work for the thermal freeze-in two-body decay. Combining with the DM relic abundance, we can then estimate the contribution to ∆N ν in light of the hidden light particle mass, but without a sensitive dependence on the decaying mass and the feeble coupling, both of which are directly responsible for the DM production. It should be mentioned that, if the antiparticleχ is not identical to χ, we should multiply Y χ,ρ(n) by a factor of two to account for the total energy (number) densities of particle and antiparticle, assuming no significant production of a χ-asymmetry. In addition, if the scalar φ is a gauge multiplet at the symmetric phase, the collision term C χ,ρ(n) should be multiplied by the number of gauge components. Nevertheless, the ratio R χ defined by Eq. (16) remains unchanged in both of these two cases.
The current CMB+BBN+Y p +D data sets a limit on ∆N ν at the level of O(0.1) [2], while the future sensitivity from BBN+CMB-S4 is expected to be ∆N ν = O(0.01), depending on the sky fraction f sky [33,34]. With these sensitivities kept in mind, we show in Fig. 1 the variation of the ratio R χ with respect to the mass of the light hidden species χ, where the magenta band is obtained by varying the freeze-in temperature from the electroweak scale (corresponding to g * ,χdec = 106.75) down to 10 MeV (corresponding to g * ,χdec = 10.76 [46]). It can be seen that, if the light species χ makes up all the DM relic density, generating a shift of ∆N ν = 0.1 requires a DM mass at m χ O(50) eV, while a shift of ∆N ν = 0.01 suggests that m χ 0.6 keV. In addition, as visualized already by Eq. (16), the shift ∆N ν would become smaller when the DM mass increases. For these sub-keV single DM candidates, they are expected to play the role of warm DM [1]. In this case, it should be mentioned that, although the mass bounds from the Lyman-α forest cannot be directly applied to the non-thermal warm DM production scenarios, these candidates are generically expected to have a lower mass bound of O(1) keV [20,[27][28][29][30][31][32]. This implies that, for a single warm DM produced via the thermal freeze-in two-body decay, the contribution to ∆N ν is negligible and hence can be probed neither by the current CMB+BBN+Y p +D nor the future BBN+CMB-S4 sensitivity, no matter what the feeble coupling strength and the decaying particle mass (m φ 1 MeV) are. Although the current and future BBN precision cannot be applied as a sensitive avenue to constrain or test the scenario where the light species χ makes up all the DM relic density, the future BBN+CMB-S4 sensitivity could be powerful in testing the mixed cold-plus-warm DM scenarios [21][22][23][24][25], where the Lyman-α constraints become weaker if the particle χ makes up only a small portion of the current DM relic density [22,23,25]. In this context, the simple scaling of the ratio R χ can serve to correlate the information of the warm DM component with the observation of ∆N ν excess in the future. As an example, we show in Fig. 1   Note that such a non-thermal warm DM component cannot be simply ruled out by the cosmological constraints discussed in Refs. [22,23,25]. The reason is that the warm DM component was presumed therein to have a thermal or a thermal-like momentum distribution, while for the non-thermal DM considered here, the momentum distribution f χ (p) is a non-trivial function of momentum and it also depends on the unknown mass spectrum m φ,ψ . As a consequence, there is no definite lower mass bound on the non-thermal sub-component DM due to its dependence on both ω χ and f χ (p). Further comprehensive numerical analyses in light of astrophysical and cosmological observables are beyond the scope of this work and will not be concerned here. Rather, we would like to highlight a simple application with respect to the ratio R χ : once a small warm portion of the DM relic density is observed, e.g., via its effect on the small-scale structures, together with a possible ∆N ν excess, we can then use the ratio R χ as an economic criterion to infer whether the origin of these two phenomena comes from a sub-keV warm DM component generated by the thermal freeze-in two-body decay.

V. VALIDITY AND APPLICATION OF Rχ
The simple scaling of the ratio R χ given by Eq. (16) is a common feature expected for the thermal freeze-in two-body decays, as well as the general situations where the coupling constant governing the DM production can be factored out from the Boltzmann integrals and the mass spectrum satisfies the hierarchy m φ m ψ,χ . If the DM production is dominated by a single channel, we can generally factor out the feeble coupling so that the ratio R χ is independent of this information, while the extraction of the decaying mass is possible only when it is the primary scale involved in the production channel. If m ψ is comparable to m φ , on the other hand, there would be an additional scale-dependence arising from the integration of the energy (number) yield Y χ,ρ(n) , and the ratio R χ would subsequently depend on the measure of the phase-space closure in the decay: Noticeably, this additional mass spectrum also appears in the first [47] and the second moment [48] of the DM momentum distribution f (p), which are defined, respectively, by Furthermore, it has been found that decreasing both of these two moments via a quasi-degenerate mass spectrum m φ ≈ m ψ can dilute the Lyman-α constraints [47,48]. Intuitively, this corresponds to transferring the energy of the decaying particle φ mostly to that of ψ so that the DM momentum inherited from the two-body decay is suppressed, making therefore the keV warm DM colder. Nevertheless, the phase-space suppression in this case also reduces the ratio R χ , as can be seen from Eq. (17). Thus, to have a significant shift of ∆N ν by a single sub-keV DM, we cannot simply evade the Lyman-α constraints with a quasi-degenerate spectrum m φ ≈ m ψ , since the contribution to ∆N ν is simultaneously suppressed. In this way, we confirm therefore the conclusion made already in Refs. [19,20] that a single warm DM candidate under the Lyman-α constraints cannot generate an observable shift of ∆N ν . It should be emphasized again that the application of the Lyman-α constraints discussed in Refs. [19,20] necessitates an explicit determination of the DM momentum distribution f (p). Certainly, the calculation of ∆N ν would be straightforward once f (p) is solved from the unintegrated kinetic equations. In the simple algorithm presented here, however, we start from the integrated Boltzmann evolution. As shown above, this allows a fast and intuitive estimate of ∆N ν in terms of the DM information itself (i.e., the DM mass and relic density), but without a sensitive dependence on the unknown production scale and coupling. Furthermore, such a freeze-in independence property could also be applied as an economic criterion to check, whether a small portion of the DM relic density from a sub-keV warm component is responsible for a possible ∆N ν excess in the future.
We would like to mention that, if the DM has a thermal history, the relic particle-number and energy densities can be easily predicted with a thermal distribution. It is then straightforward to obtain that the ratio R χ in this case has a similar form to Eq. (16) up to a small numeric correction. For generating an observable shift of ∆N ν from a light thermal DM, the Lyman-α constraints will directly come into play. Consequently, the single thermal warm DM is readily ruled out, while the thermal component in the cold-plus-warm DM scenarios would face a strong constraint from the Lyman-α forest data.
As a final remark, we should emphasize that the freezein independence is a result of the thermal freeze-in mechanism, where the scalar φ can decay into and scatter with the SM particles in realizing the φ thermalization in the early Universe. Given a great deal of DM scenarios, it is not true that the simple scaling R χ specified by Eq. (16) always exists. For instance, if φ is also a component of the hidden sector where some dark Z 2 symmetry resides at the φ − χ sector, the super-WIMP mechanism [49] can be at work to generate χ in late-time decay after φ freezes out from the SM plasma. In this case, φ is a WIMP while χ the super-WIMP DM, and one will then get the following scaling: (20) where Y χ,n is the relic abundance inherited from the yield Y φ,n after φ has decayed away to χ. Y φ,n is determined at the frozen-out temperature T φ,f o , which depends on the particle information of φ, including its gauge representation, interactions and mass. y 2 is extracted from the amplitude squared in Eq. (11), while F is a calculable function in which f φ (y, m φ ) is no longer a thermal distribution function due to the late-time out-of-equilibrium decay and hence depends on the decay coupling y and the decaying scalar mass m φ . This indicates that the ratio Y χ,ρ /Y χ,n will depend on (y, m φ ) in a non-trivial way. Therefore, if φ and χ form a hidden dark sector such that the super-WIMP mechanism plays the role for χ production, the simple scaling of R χ in terms of the DM mass cannot be realized.

VI. CONCLUSION
If the hidden light particle has a relativistic history at the BBN epoch while becomes non-relativistic prior to the CMB regime, it is in principle possible to probe these light species via the BBN precision measurements. Motivated by this, we have demonstrated in this work that, in the thermal freeze-in production regime, the current and future BBN sensitivities cannot probe a single keV-scale warm DM under the severe constraint from Lyman-α for-est data, since the latter renders only a negligible contribution to ∆N ν . However, the future forecast from CMB-S4, when combined with the BBN measurements, has the potential to test the warm component of cold-plus-warm DM scenarios. Such a simple observation presented here depends only on the particle-number and the energy evolution, the ratio of which can be applied to predict the shift ∆N ν , basically without the sensitive information of the DM freeze-in scale and feeble coupling.
While the simple scaling of R χ given by Eq. (16) exists in the standard freeze-in mechanism, it is not always true for some other interesting DM scenarios. Nevertheless, we expect that this project deserves further studies in the future, as it brings us an interesting and simple avenue to probe the light hidden species that are otherwise difficult to probe directly at terrestrial colliders and indirectly in flavor physics.