Observational Constraints on Secret Neutrino Interactions from Big Bang Nucleosynthesis

We investigate possible interactions between neutrinos and massive scalar bosons via $g^{}_{\phi} \overline{\nu} \nu \phi$ (or massive vector bosons via $g^{}_V \overline{\nu} \gamma^\mu \nu V^{}_\mu$) and explore the allowed parameter space of the coupling constant $g^{}_{\phi}$ (or $g^{}_V$) and the scalar (or vector) boson mass $m^{}_\phi$ (or $m^{}_V$) by requiring that these secret neutrino interactions (SNIs) should not spoil the success of Big Bang nucleosynthesis (BBN). Incorporating the SNIs into the evolution of the early Universe in the BBN era, we numerically solve the Boltzmann equations and compare the predictions for the abundances of light elements with observations. It turns out that the constraint on $g^{}_{\phi}$ and $m^{}_\phi$ in the scalar-boson case is rather weak, due to a small number of degrees of freedom. However, in the vector-boson case, the most stringent bound on the coupling $g^{}_V \lesssim 6\times 10^{-10}$ at $95~\%$ confidence level is obtained for $m^{}_V \simeq 1~{\rm MeV}$, while the bound becomes much weaker $g^{}_V \lesssim 8\times 10^{-6}$ for smaller masses $m^{}_V \lesssim 10^{-4}~{\rm MeV}$. Moreover, we discuss in some detail how the SNIs affect the cosmological evolution and the abundances of the lightest elements.


Introduction
As one of the pillars of the standard model of cosmology, Big Bang nucleosynthesis (BBN) opens a unique window to the early Universe [1][2][3][4][5][6] and new physics beyond the standard model of elementary particles [7][8][9]. Based on the standard models of both cosmology and particle physics, the theory of BBN itself essentially contains only one free parameter, i.e., the baryon-to-photon density ratio, which has been determined to be η ≡ n b /n γ = (6.047 ± 0.074) × 10 −10 from the precision measurement of the cosmic microwave background (CMB) by the Planck collaboration [10]. The observed ratio between the primordial abundance of deuterium and that of hydrogen D/H| p = (2.53 ± 0.04) × 10 −5 [11][12][13][14], together with the primordial mass fraction of 4 He, i.e., Y p ≡ ρ( 4 He)/ρ b = 0.2449 ± 0.0040, indicates that 5.7 × 10 −10 η 6.7 × 10 −10 at 95 % confidence level (CL) [15,16], which is remarkably consistent with the CMB determination of the baryon density. § Therefore, any new physics leading to significant deviations from the standard BBN predictions for the light element abundances will receive restrictive constraints.
In this work, we investigate the observational constraints from BBN on the secret neutrino interaction (SNI) with a massive scalar or vector boson. More explicitly, we consider the SNI only for the left-handed neutrino fields and the relevant Lagrangian can be written as where φ and V µ are the fields for the scalar and vector boson with masses m φ and m V , respectively. For simplicity, the coupling constants g αβ φ and g αβ V are assumed to be both flavor diagonal and universal for three neutrino flavors, namely, g αβ φ = g φ δ αβ and g αβ V = g V δ αβ . The main motivation for such an investigation is two-fold. First, the interaction among neutrinos themselves has never been directly tested in terrestrial experiments, since neutrinos participate only in the neutralcurrent weak interaction in the standard electroweak theory and there has not been an attempt to collide two neutrino beams. In contrast, the early Universe and the core-collapse supernovae, where the neutrino number density is extremely high and the interaction among themselves is important, serve as ideal places to constrain the SNI [18][19][20][21][22]. Second, the SNI is expected in many particle-physics models, which have been proposed to generate tiny neutrino masses [23][24][25][26][27][28] or solve the potential problems associated with dark matter [29].
Stringent constraints on the SNI with a massless or massive scalar boson have been derived from observations of the CMB and cosmological large-scale structure formation [30][31][32][33][34][35][36][37][38][39][40], the supernova SN1987A [41][42][43][44][45], and other experiments [46][47][48][49][50][51]. These experimental constraints are also applicable with some modifications to the case of a vector boson. To this end, a detailed study of the BBN bounds on the SNI in both the scalar and vector cases is lacking, except for a brief discussion in Ref. [52]. For this reason, we now extend the previous work by incorporating the SNI into the cosmological evolution during the BBN era and examining its impact on the light element abundances. § The BBN theory can also predict the primordial abundances of 3 He and 7 Li. However, the only data on 3 He are available for the solar system and the high-metallicity regions in our Galaxy and it is difficult to infer its primordial abundance [15]. On the other hand, the observed relative abundance of lithium is 7 Li/H| p = (1.6 ± 0.3) × 10 −10 , showing a discrepancy in the baryon density between the BBN and CMB estimates [9,17]. Since the lithium abundance remains an unresolved issue, it will not be used to draw any constraints in this work. Figure 1: The Feynman diagrams for the decay φ → ν + ν (a), the annihilation ν + ν → φ + φ (b1) and (b2), the elastic scattering ν + φ → ν + φ (c1) and (c2) processes in the presence of SNI with a massive scalar boson φ. The corresponding diagrams in the vector case can be obtained by replacing φ with V accordingly.
The remaining part of this work is organized as follows. In Sec. 2, we set up the general theoretical framework to study the SNI in the BBN era. A general discussion on how the presence of new particles and interactions affects the standard BBN theory is given. Then, in Sec. 3, we numerically solve the Boltzmann equations for the cosmological evolution and the nucleosynthesis of light elements, where the BBN constraints on the coupling constant and the mass are derived. Finally, in Sec. 4, we summarize our main results and draw our conclusions.

Simple Arguments
As mentioned below Eq. (1), the coupling constant g φ or g V between neutrinos and the new particle φ or V is assumed to be flavor conserving and universal. The relaxation of such an assumption may lead to slight differences. For instance, if φ or V is coupled exclusively to ν e and copiously produced after the decoupling of ν µ and ν τ , the energy density of the decoupled ν µ and ν τ will not be modified, rendering the constraint on the coupling constant relatively weaker. For simplicity, we ignore the flavor dependence of the SNI and treat neutrino flavors equally in the evolution of our Universe. ¶ Through the interaction given in Eq. (1), the scalar boson φ ¶ If neutrino flavor conversion and non-instantaneous decoupling of neutrinos are taken into account, the effective number N eff of neutrino species (at the temperature far blow 1 MeV) defined via the neutrino-to-photon ratio of energy densities ρ ν /ρ γ = N eff · 7/8 · (4/11) 4/3 will be shifted from N eff = 3 to N eff = 3.045 [53][54][55][56]. Such a tiny difference will not be important compared to the radical changes due to the new physics under consideration.
can be generated by the inverse decay ν + ν → φ and the neutrino-antineutrino annihilation ν + ν → φ + φ, for which the Feynman diagrams are shown in Fig. 1. For the vector boson V , we have basically the same production processes. As pointed out in Ref. [52], when φ or V becomes in thermal equilibrium around the temperature T 1 MeV, it contributes to the extra radiation represented in terms of ∆N eff ≡ N eff − 3, where the effective number of neutrino species is defined as N eff ≡ (ρ r − ρ γ )/(ρ std ν /3) with ρ r and ρ std ν being the energy density of all radiation and the neutrino energy density in terms of the photon temperature T γ in the limit of instantaneous decoupling, respectively. With this definition, N eff in the standard case without SNI will be fixed to three during the whole BBN era. In the scalar case with m φ 1 MeV, we have only one extra relativistic degree of freedom, corresponding to ∆N eff = 1/2 · 8/7 ≈ 0.57, whereas in the vector case with m V 1 MeV, we have three helical states, indicating ∆N eff = 3/2 · 8/7 ≈ 1.71. It is worthwhile to notice that the temperatures of photons and neutrinos remain equal before the electron-positron annihilation at T 0.511 MeV. If m φ or m V is much larger than 1 MeV, φ or V is non-relativistic and its number density will be suppressed by the Boltzmann factor e −m φ /T or e −m V /T , significantly reducing the contribution to ∆N eff . Then, we require that the upper bound ∆N eff < 1 at 95 % CL [57] should be satisfied for T 1 MeV. Obviously, there are essentially no constraints on g φ and m φ , whereas the constraints on g V and m V could be restrictive.
In the above discussion, we have only considered the situation in which V can be in thermal equilibrium at T 1 MeV, but whether this is the case or not depends on g V and m V . On the other hand, even if V cannot be brought into thermal equilibrium, it still contributes to the total energy density of our Universe. Therefore, the Boltzmann equations for the distribution functions of both neutrinos and V should be implemented in the latter case to calculate ∆N eff . Before going into details of the Boltzmann equations, we make some comments on the constraints on g V and m V by simple dimensional analysis: • Around T 1 MeV, the Universe is certainly radiation-dominated, so the Hubble expansion rate is given by H ≈ 1.66 √ g * T 2 /M Pl , where g * = 10.75 denotes the effective number of relativistic degrees of freedom and M Pl 1.221 × 10 19 GeV is the Planck mass. For V to be thermalized, we have to calculate its production rate Γ V and demand Γ V H at T 1 MeV.
• For a relatively large mass m V 1 MeV, the inverse decay ν + ν → V could be quite efficient, since the decay rate in the rest frame is proportional to both g 2 V and m V . As an order-of-magnitude estimate, we obtain Γ D where the last factor is the Lorentz factor E /m V ≈ 3T /m V , arising from a boost to the comoving frame. When the inverse decay dominates the production and brings V into thermal equilibrium, i.e., Γ D V H, we arrive at Consequently, if V is thermalized at T 1 MeV via the inverse decay, the upper bound ∆N eff < 1 can be translated into g V 2.2 × 10 −10 (1 MeV/m V ). This is well consistent with the result from a more detailed calculation in Ref. [52].
• For an extremely small mass m V 1 MeV, the inverse decay becomes inefficient and the annihilation ν + ν → V + V will take over in thermalizing V . Since the masses of neutrinos and V are sufficiently small, the only relevant energy scale in question is just T . Hence, Γ V can be estimated to Similarly, the upper bound ∆N eff < 1 will restrict g V into the region of g V 4.6 × 10 −6 . An accurate calculation of the total energy density or ∆N eff in the case of m V 1 MeV has not yet been performed in the literature.
The exact calculation of the total energy density for an arbitrary coupling constant g φ (or g V ) and an arbitrary mass m φ (or m V ) calls for the implementation of Boltzmann equations. First, both neutrinos and the new boson φ or V could deviate from the thermal distributions, so the determination of ∆N eff requires numerical solutions to the true distribution functions. Second, since the BBN takes place for a wide range of temperature (e.g., from T = 1 MeV to T = 0.01 MeV), a naive requirement for ∆N eff < 1 at T 1 MeV oversimplifies the picture of relevant physics.

Boltzmann Equations
In order to fully take into account the decay, annihilation, and scattering processes shown in Fig. 1, we need to solve a complete set of Boltzmann equations for the distribution functions of neutrinos and the new particle φ or V . The general theoretical framework for the cosmological evolution can be found in a number of excellent books on cosmology (see, e.g., Refs. [58][59][60]). Therefore, we only outline the strategy for our computations.
First, the Hubble parameter H(t) ≡ȧ(t)/a(t) is governed by the Friedmann equation , where a(t) is the scale factor and ρ is the total energy density. The evolution of the energy density satisfiesρ(t) = −3H(ρ + P ), where both ρ and pressure P can be solved for the given particle contents and their distributions. As ρ = ρ γ + ρ ν + ρ e + ρ b + ρ φ/V is still dominated by radiation in the BBN era, we can safely ignore the contribution ρ b from baryons and count all those from photons (γ), neutrinos (ν), electrons (e), and the new boson φ (or V ). This sets up the evolution of the cosmological background.
Second, in the homogeneous and isotropic Universe with the Friedmann-Lemaître-Robertson-Walker metric, the relevant distribution functions f i (|p i |, t) for i = ν, ν, and φ (or V ) fulfill the following Boltzmann equations [61]: where the quantities C i D , C i A , and C i E stand for the collision terms of the decay, annihilation, and elastic scattering processes, respectively. In fact, the last term C i SM in the right-hand side of Eq. (4) collectively includes all the relevant scattering processes for neutrinos in the standard model, such as νν ↔ νν, νν → νν, νν ↔ e + e − , and νe − → νe − , where the neutrino (antineutrino) flavor indices have been suppressed. Since the neutrino interactions in the standard model have been extensively studied in the literature [62-64], we will not explicitly show them in Eq. (4), but have indeed included them in our numerical calculations. For the SNI part, assuming the scalar boson φ, we have where dp is the Dirac delta function for four-momentum conservation. The matrix elements squared |M D | 2 , |M A | 2 , and |M E | 2 are averaged over the initial and final spins. For the annihilation process, one should take care of the symmetric factors, arising from the identical particles in the initial and final states, and the changes of particle numbers in each specific process. Furthermore, the contribution from the elastic scattering between ν and φ is not explicitly shown in C φ E , but should be added. In fact, for the Boltzmann equations of neutrinos and antineutrinos, we have also included the standard model processes, which establish a thermal contact between the neutrino sector and the system of photons, electrons, and baryons.
Finally, we come to the nuclear reactions for the generation of light elements. At high temperatures, both neutrons (n) and protons (p) are in thermal equilibrium, so the neutron-to-proton ratio is given by n/p = e −(m n −m p )/T . After the weak interaction freezes out and neutrinos decouple from the thermal bath at T 1 MeV, we have n/p ≈ 1/6, which will drop to 1/7 by the time of nuclear reactions due to beta decays of free neutrons. The first process is the formation of deuterium (D) via p + n → D + γ, which is at work efficiently after the photo-disintegration rate is suppressed at T 0.1 MeV. As neutrons will ultimately be integrated into the most stable light element 4 He, one can estimate its mass fraction via Y p = 2(n/p)/(1 + n/p) ≈ 0.25 [15]. Although Y p is not quite sensitive to the nuclear reactions, the abundances of D, 3 He, and 7 Li relative to that of H are of the order of 10 −5 or even smaller and will be greatly affected by the detailed reactions. In order to numerically calculate Y p and D/H| p , we implement the publicly available code AlterBBN of Ref. [65], which is actually based on the original Fortran code first presented in Refs. [66,67] and updated with the latest cross sections of relevant nuclear reactions. Another well-known code PArthENoPE has also been widely used [68,69]. Note that the neutron lifetime τ n = (880.3 ± 1.1) s will be used in our calculation as an input value [16].

Extra Radiation
As a numerical support for the simple arguments given in Sec. 2.1, we now attempt to compute the total energy density by solving the Boltzmann equations for the distribution functions of neutrinos and the new particle φ or V and to determine the extra radiation in terms of ∆N eff at T 1 MeV. For the moment, the nucleosynthesis of light elements is not initiated. We start by specifying the initial conditions at a high temperature T = T γ = T ν = 10 MeV. For later convenience, a(t) is normalized to 1/T γ and two dimensionless parameters x ≡ ma and q ≡ |p|a are introduced, where m can be an arbitrary mass scale and is set to 1 MeV in practice. Photons, neutrinos, and electrons initially follow the distributions in thermal equilibrium with zero chemical potentials, while the initial abundance of φ or V is assumed to be negligible. With the help of the two dimensionless parameters, we can recast the Boltzmann equations into the following form and obtain xdρ/dx = −3(ρ + P ), where ρ is in general composed of two parts: the thermal-bath sector ρ γ and ρ e and the neutrino sector ρ ν and ρ φ (or ρ V ). The former sector can be described by the equilibrium distribution with T γ . For the latter sector, the energy density should be calculated from the real distribution functions f ν and f φ (or f V ). There is essentially no difference between neutrinos and antineutrinos, so ρ ν actually represents the energy density of both. Then, it is straightforward to derive Solving Eqs. (8) and (9) with H = 8πρ/(3M 2 Pl ), we obtain the total energy density of radiation for any given values of g φ (or g V ) and m φ (or m V ) and can extract the extra radiation ∆N eff .
In Fig. 2, the numerical results are presented, where the contours of ∆N eff have been plotted in the plane of m φ and g φ in the scalar-boson case (left panel) and in the plane of m V and g V in the vector-boson case (right panel). Compared with the previous study in Ref. [52], the parameter space has now been extended to the region of much smaller masses, for which the annihilation processes become dominant in the production of φ or V . Some interesting features can be observed from Fig. 2. First, in both plots, one can see that the contours turn out to be flat at the small-mass end, e.g., m φ or m V ≈ 10 −5 MeV, indicating that the results of ∆N eff can be simply extrapolated to the cases of even smaller masses. Second, at the high-mass end, even if φ or V can be in thermal equilibrium, the Boltzmann factor leads to a suppression of ∆N eff . This is why ∆N eff decreases quickly when the mass increases. Third, in between low and high masses, the product of g φ m φ or g V m V is nearly constant for a given ∆N eff , which can be understood by the simple estimate in Eq. (2). However, the maximum of ∆N eff in the scalar case is 0.57, when φ is relativistic and the SNI can bring it into thermal equilibrium with neutrinos. As a consequence, the BBN bound ∆N eff 1 from Ref. [57] does not have any constraining power on g φ and m φ .

Light Element Abundances
After having numerically computed the extra radiation, we proceed with the real evolution of light element abundances by combining the Boltzmann equations with the AlterBBN code [65]. Solving the Boltzmann equations numerically, we can obtain the necessary information on the evolution of the cosmic background. Such an information will be input to AlterBBN as an alternative model of cosmology, and thus, the light element abundances can be calculated at any instant of time or temperature.
As have already been observed in previous sections, φ cannot contribute significantly to an extra radiation during the BBN era. In order to illustrate the main effects of the SNI in the BBN, we will concentrate on V and try to capture the most important points by analyzing four different cases of couplings and masses.
Case I: g V = 10 −4 and m V = 10 −5 MeV. In Fig. 3 (a1), we show the evolution of the comoving energy densities of γ's, ν's, and V 's. In Fig. 3 (a2), the neutron-to-proton ratio n/p and the mass fraction of helium Y p are given. The standard theory (dashed curves) is rather simple, which means that the comoving energy density of γ's increases due to the heating from electron-positron annihilation, while that of ν's remains nearly unchanged, since ν's are almost decoupled from γ's after 1 MeV. In contrast, in this non-standard case, the coupling constant g V = 10 −4 is so strong that the V 's have been tightly coupled with ν's even before T γ = 10 MeV or a = 0.1/MeV. The non-standard (solid blue curve) comoving energy density of γ's is smaller than the standard one (dashed blue curve), since the entropy has been transferred into V 's that are abundantly produced. The same behavior of the comoving energy density of ν's (green curves) is also found. However, one should notice that the physical processes actually depend on T γ or ρ γ instead of a(t). One can also observe that the energy density ratio of the extra degree of freedom to γ's by far exceeds that of the standard case. Namely, the extra radiation in terms of the effective number of neutrino species is N eff = 4.71. Consequently, this would accelerate the expansion of the Universe, leading to an earlier decoupling of the weak interaction that interchanges neutrons and protons and a larger value of n/p would then remain for 4 He synthesis. This can been observed in Fig. 3 (a2), where the asymptotic value of Y p is 8.5 % larger than the standard value. For comparison, the 1σ error of the observed Y p is only 1.6 %. Therefore, this case should give a restrictive constraint.
Case II: g V = 10 −7 and m V = 0.05 MeV. In Fig. 3 (b1) and (b2), the numerical results for  Figure 3: The cosmological evolution of the comoving energy densities for neutrinos ν (green curves), photons γ (blue curves), and vector bosons V (red curves) is presented in the left panel and that of the neutron-to-proton ratio n/p (brown curves) and the helium mass fraction Y p (purple curves) in the right panel. The plots (a1) and (a2) are given for g V = 10 −4 and m V = 10 −5 MeV in Case I, while (b1) and (b2) for g V = 10 −7 and m V = 0.05 MeV in Case II. Note that the evolution of the comoving energy densities is plotted with respect to the scale factor a(t) that has been normalized to 1/T γ at very high temperatures. The dotted vertical lines at a = 1/MeV or T γ = 1 MeV in the plots (a1) and (b1) represent basically the epoch of neutrino decoupling. The dashed curves are for the standard theory, while the solid curves for the non-standard one. The baryon-to-photon ratio η takes on the value, which minimizes the χ 2 function for each case.
this case are given. Since g V = 10 −7 is now much smaller compared to that in Case I, the comoving energy densities of γ's (blue curves) and ν's (green curves) coincide with those of the standard values at T γ = 10 MeV. Shortly later on, one can observe sizable deviations due to the production of V 's. Around a ≈ 50/MeV or T γ ≈ 0.02 MeV when the V 's remain in thermal equilibrium with ν's, the comoving energy density ρ V (red curves) is suppressed by the Boltzmann factor, since m V = 0.05 MeV > T γ ≈ 0.02 MeV. The reduction of the number density of V 's is taking place via both the decay V → ν + ν and the annihilation V + V → ν + ν processes. This is the reason Figure 4: The cosmological evolution of the comoving energy densities for neutrinos ν (green curves), photons γ (blue curves), and vector bosons V (red curves) is presented in the left panel and that of the neutron-to-proton ratio n/p (brown curves) and the helium mass fraction Y p (purple curves) in the right panel. The plots (a1) and (a2) are given for g V = 2 × 10 −7 and m V = 0.003 MeV in Case III, while (b1) and (b2) for g V = 10 −6 and m V = 2.5 MeV in Case IV. The notations and other input parameters are the same as for Fig. 3.
why the comoving energy density of ν's takes over that of γ's after these processes are almost completed. Since the V 's are thermally populated already by the time of neutrino decoupling at T γ = 1 MeV, ∆N eff is the same as in Case I and its impact on n/p and Y p is also quite similar.
Case III: g V = 2 × 10 −7 and m V = 0.003 MeV. In Fig. 4 (a1) and (a2), the numerical results for this case are displayed. From Fig. 4 (a1), one can observe that the V 's (red curves) are mostly generated after the neutrino decoupling at T γ ≈ 1 MeV, rendering the energy density of γ's (blue curves) to be nearly unchanged. In this case, the extra radiation is found to be ∆N eff ≈ 0.5. However, T ν is severely reduced by thermalizing the V 's. This picture cannot be effectively described by ∆N eff . As indicated in Fig. 4 (a2), the increase of n/p and Y p is more significant than those in the two previous cases, although ∆N eff is smaller. The true reason is that the lower T ν would cause an earlier decoupling of the weak interaction, implying a larger value of n/p. In this case, we find that the final Y p deviates from the standard value by 10 %. This actually results in the peak of ∆χ 2 on the edge of the excluded region in Fig. 6, which will be discussed in the next subsection.
Case IV: g V = 10 −6 and m V = 2.5 MeV. In Fig. 4 (b1) and (b2), the numerical results for this case are shown. Due to the large mass, the V 's are mainly produced via the inverse decay ν + ν → V and have been tightly coupled to ν's before T γ = 10 MeV. In this case, the extra radiation is about ∆N eff ≈ 1.5. However, the light element abundances seem to be unaffected, as indicated in Fig. 4 (b2). This can be interpreted as a cancellation between the effects of a higher T ν and a larger expansion rate. One may extend the discussion on this case to future works, since this might provide a possible solution to the 7 Li abundance problem [9], which is, however, beyond the scope of the present work.
In Fig. 5, the evolutions of g * and N eff as functions of T γ have been displayed for different parameters in the vector-boson case. As we have mentioned before, N eff has directly been calculated from the true distribution functions of ν's and V 's, which themselves are found by solving the relevant Boltzmann equations. Since N eff itself does evolve with respect to temperature, it is not correct to just draw the constraints on g V and m V by simply putting an upper bound on N eff at any instant of temperature or time. Although the V 's contribute to N eff by at most 1.71 when they are relativistic and in thermal equilibrium, N eff ≈ 5 can be reached in Cases II and IV, as shown in the right panel of Fig. 5. The main reason is that the V 's decay into ν's, when T γ drops below m V , and transfer their energies into ν's. This feature can also be observed from Figs. 3 (b1) and 4 (b1). Since the V 's are non-relativistic at this stage, their energy densities will ∆χ 2 = 1.8 be diluted less significantly than those of relativistic particles, such as ν's. Therefore, the total energy density in terms of N eff could go slightly beyond the value 4.71. For comparison, in the left panel of Fig. 5, we also present g * (T γ ), which is defined through the radiation energy density ρ r (T γ ) = π 2 g * (T γ )T 4 γ /30. From the above discussion of the four different cases, it is now clear how the light element abundances are affected by the SNI with a massive scalar or vector boson. The impact of a scalar boson is moderate, since it has a smaller number of degrees of freedom compared to a vector boson. In general, the introduction of the SNI will produce additional radiation to accelerate the expansion of the Universe. In addition, the thermal contact between the new particle φ or V and neutrinos will change T ν after neutrino decoupling. Both effects will be important for the evolution of the light element abundances. It should be noticed that the abundance of deuterium is also modified when the SNI is present, but the deviations from the standard theory for Cases I-IV are within 1 % when the χ 2 function is minimized with respect to η. Thus, the evolution of deuterium is not included. Details on the minimization of χ 2 can be found in the next subsection, where more discussions on the evolution of light element abundances will be given.

Final Constraints
The ultimate goal of this work is to constrain the parameter space by the observations of the primordial abundances of deuterium and helium. For this purpose, we have to confront the theoretical predictions of these quantities with the observed values. With the help of our modified version of the AlterBBN code, we are able to calculate the light element abundances Y i (η, g s , m s ) given the baryon-to-photon ratio η, the coupling constant g s , and the mass m s . Here the subscript "s" denotes either φ for the scalar boson or V for the vector boson. Apart from these input parameters, we need to take into account the theoretical errors on the nuclear reaction rates and the neutron lifetime. This is important, since the observational errors are currently comparable to the theoretical ones. To this end, we adopt the simple method of linear error propagation from Ref. [70] and define the χ 2 function as where S ij ≡ (σ 2 th ) ij + (σ 2 ex ) ij is the covariance matrix of squared errors and the indices i and j refer to the light elements of our interest. We only consider the well-measured abundances of deuterium and helium, i.e., i and j run over D and 4 He. In Eq. (10), Y th i and Y ex i are the theoretical predictions and the experimental values of the abundances, respectively. In our calculations, the theoretical errors (σ 2 th ) ij are estimated by using the method of linear error propagation. For simplicity, we just consider the errors of the twelve most relevant nuclear reactions [70]. The theoretical errors on those nuclear reactions can be found in Refs. [65,68,71], where the astrophysical S factors and their corresponding uncertainties are quantified as polynomial functions of temperature. One should be referred to Refs. [72][73][74][75][76][77][78][79][80][81][82] for recent developments in this aspect. On the other hand, the experimental errors are assumed to be uncorrelated, namely, (σ 2 ex ) ij = δ ij σ i ex σ j ex , where the standard deviation σ i ex is assumed for the corresponding individual observable. In Ref. [11], the newly added He I λ10830 infrared emission lines have been incorporated together with the traditional visible lines from the metal-poor ionized hydrogen regions of compact blue galaxies to obtain Y p = 0.2449 ± 0.0040. The relative abundance of deuterium D/H| p = (2.53 ± 0.04) × 10 −5 is inferred from the observation of high-redshift interstellar clouds in Ref. [12], in which all existing data have been systematically studied.
Then, we are ready to compute χ 2 (η, g s , m s ) by scanning over the parameter space of η, g s , and m s . Since we are interested in the BBN constraints on g s and m s , the values of χ 2 (η, g s , m s ) are minimized with respect to η in the range from 10 −9 to 10 −10 . Since the SNI may also affect the formation of the CMB, the observational information on η from the CMB data is not used for self-consistency. Since the abundance of deuterium is very sensitive to η and the one of 4 He is not, minimizing χ 2 with respect to η is equivalent to minimizing the deuterium deviation with respect to η. This is also why the deviations of the abundance of deuterium from the standard values are quite small in Fig. 3 and Fig. 4 when the minimization is performed. In Fig. 6, the final results are presented, where the contours and the density map of the χ 2 function have been plotted for the scalar-boson case (left panel) and the vector-boson case (right panel). In general, the excluded regions are similar to those in Fig. 2, which are basically described by simple arguments of ∆N eff . Nevertheless, some special features in Fig. 6 can only be explained after a full calculation of light element abundances is accomplished.
In the left panel of Fig. 6, there is a gap between the small scattered circles and the boundary of a large connected region, corresponding to the contours of ∆χ 2 = 1.8. The reason is essentially the same as already mentioned in Case III of the previous subsection. It is the freezing-in of φ before or after the neutrino decoupling that makes the key difference. In the former case, when φ is in thermal equilibrium before the neutrino decoupling, a large value of g φ is required. Thus, T γ and T ν must be changed simultaneously such that their ratio is almost the same as that of the standard case. Hence, the light abundances are affected by the extra radiation, or equivalently, a large expansion rate. However, in the latter case, when φ freezes in after the neutrino decoupling, a smaller g φ is needed. Although the contribution to the extra radiation is small, the production of φ's from ν's will reduce T ν , leading to an earlier decoupling of the weak interaction. Therefore, for a fixed m φ , both a large coupling and a small one can give rise to the same χ 2 function, successfully explaining the observed features. In the right panel of Fig. 6, the gap does not show up, but one can notice a dark region along the 99 % contour in the density map, which has the same origin as the appearance of the smaller circles in the left panel.
With the above detailed calculations, we now make some remarks on the flavor dependence of the SNI and discuss its impact on the final observational constraints. If the mediator φ or V is only coupled to muon and tau neutrinos (antineutrinos), the influence on BBN will be just the modification of the total energy density, which can entirely be represented by N eff . However, if the mediator is coupled to the electron neutrinos (antineutrinos), the average energy or temperature of ν e (ν e ) will be changed via the production of the mediator, which will affect the total energy density and directly modify the rates of weak interactions ν e +n ↔ p+e − and ν e +p ↔ n+e + that are responsible for the n/p ratio. As for the final constraints on the couplings and the mediator masses, the main difference between the electron and muon (or tau) neutrinos can be clearly observed by comparing the numerical results shown in Fig. 6 and those in Fig. 2. In the former case, the dark regions along the indicated contours in Fig. 6 signify the additional impact on the n/p ratio in the case of the SNI for ν e .

Summary and Conclusions
We have performed a detailed study of the SNI in the BBN era. In order to better understand the production and the evolution of the involved new scalar or vector boson, we have implemented the Boltzmann equations for the distribution functions of neutrinos and this new particle. Furthermore, these Boltzmann equations have been combined with the AlterBBN code to analyze the effects of the SNI in the evolution of light element abundances. The observed primordial abundances of deuterium and helium have been used to constrain the coupling and the mass of the new particle. Such a study is very helpful in exploring the intrinsic properties of massive neutrinos, including the origin of neutrino masses and the interactions among neutrinos, and in constraining other new-physics scenarios beyond the standard model.
In the scalar-boson case, it has been demonstrated that the BBN bound on g φ and m φ is very weak. However, in the vector-boson case, very stringent bounds on g V can be obtained for a wide range of masses, i.e., 10 −5 MeV m V 5 MeV. The most stringent bound g V 6 × 10 −10 at 95 % CL is achieved for m V 1 MeV. The bound on g V will be significantly relaxed for much smaller masses, namely, g V 8 × 10 −6 for m V 10 −4 MeV at the same CL. The BBN constraint is comparable to the supernova bound [45], but weaker than the bounds drawn from decays of kaon and weak gauge bosons in the low-mass region [51].
It is worthwhile to emphasize that the CMB data will be very useful in constraining the SNI. For instance, although the scalar boson φ can at most contribute an extra radiation of ∆N eff = 0.57, the CMB bound on the effective number of neutrino species is N eff = 3.14 +0. 44 −0.43 at 95 % CL [10]. In addition, when the temperature drops below m φ , the Boltzmann suppression processes conserve the entropy, but not the comoving energy density, enhancing the extra radiation to ∆N eff ≈ 0.75 in the CMB epoch. Similar arguments can be applied to the vector boson V as well. However, it should be noted that the SNI in the CMB epoch may alter the power spectrum of anisotropy via a tight coupling [30][31][32][33][34][35][36][37]40]. Fortunately, both a larger N eff and a stronger SNI enhance the power spectrum of anisotropy, leading to more restrictive constraints.