SSolving the Hubble tension without spoiling Big Bang Nucleosynthesis
Guo-yuan Huang ∗ and Werner Rodejohann † Max-Planck-Institut f¨ur Kernphysik, Postfach 103980, D-69029 Heidelberg, Germany (Dated: February 9, 2021)The Hubble parameter inferred from cosmic microwave background observations is consistently lowerthan that from local measurements, which could hint towards new physics. Solutions to the Hubbletension typically require a sizable amount of extra radiation ∆ N eﬀ during recombination. However,the amount of ∆ N eﬀ in the early Universe is unavoidably constrained by Big Bang Nucleosynthesis(BBN), which causes problems for such solutions. We present a possibility to evade this problemby introducing neutrino self-interactions via a simple Majoron-like coupling. The scalar is slightlyheavier than 1 MeV and allowed to be fully thermalized throughout the BBN era. The rise ofneutrino temperature due to the entropy transfer via φ → νν reactions compensates the eﬀect ofa large ∆ N eﬀ on BBN. Values of ∆ N eﬀ as large as 0 . Introduction . —The Hubble parameter inferred fromthe Planck observations of the cosmic microwave back-ground (CMB), H = 67 . ± . / s / Mpc , is in ten-sion with that of local measurements at low red-shifts.To be speciﬁc, the result from the Hubble Space Tele-scope (HST) by observing the Milky Way Cepheids is H = 74 . ± .
42 km / s / Mpc , which exceeds the valueof Planck experiment by a 4 . σ signiﬁcance. Combiningthe HST result with a later independent determination yields a lower value of H = 72 . ± .
19 km / s / Mpc, butthe tension still persists at 3 . σ level.The tension for the Hubble parameter could suggestthe existence of new physics beyond the Standard Modelor beyond the ΛCDM framework . There is a strongpositive correlation between H and an extra radiation,∆ N eﬀ = N eﬀ − . N eﬀ during recombination one can lift the Hub-ble parameter. However, increasing N eﬀ delays matter-radiation equality and modiﬁes the CMB power spec-trum. This, in turn, can be compensated by introducingnon-standard neutrino self-interactions during recombi-nation [5, 6]. Thus, a successful particle physics modelto explain the Hubble tension needs to provide a signif-icant amount of ∆ N eﬀ in the early Universe, as well as“secret” neutrino interactions. In the original ﬁts withPlanck 2015 data, two modes with self-interacting neutri-nos of the form G eﬀ νννν are identiﬁed, which are givenin Table I. The mode with strongly interacting neutrinos(SI ν ) is basically excluded by various terrestrial experi- Parameter log ( G eﬀ · MeV ) ∆ N eﬀ η SI ν − . +0 . − . . ± .
29 6 . +0 . − . MI ν − . +1 . − . . ± .
28 6 . ± . σ ranges of two modes in theﬁt of Planck 2015 data . SI ν (MI ν ) stands for the strongly(moderately) interacting neutrino mode; η ≡ η b × rep-resents the baryon-to-photon ratio. These results are updatedwith the Planck 2018 data in Refs. [7–10]. ments [11–15], so we shall conﬁne ourselves to moderatelyinteracting neutrinos (MI ν ).Among many attempts [16–24], one of the simplestpossibilities is the Majoron-like interaction [21–24] L ⊃ g αβ φν α L ν c β L , (1)where α and β run over ﬂavors e , µ and τ , and g αβ areﬂavor-dependent coupling constants. After the neutrinotemperature drops below m φ , the interactions amongneutrinos will be reduced to an eﬀective coupling G eﬀ = | g | /m φ . The ﬂavor-speciﬁc couplings g ee and g µµ areseverely constrained by laboratory searches [11–15], andwe are only left with g ττ to accommodate the MI ν modeduring recombination.In our model, the scalar particle φ with a mass m φ increases N eﬀ . As long as the coupling in Eq. (1) isstrong enough, φ will be in thermal equilibrium withthe neutrino plasma, contributing to extra radiation by∆ N eﬀ = 1 / · / (cid:39) .
57 for m φ (cid:28) T ν , where T ν is theplasma temperature. Note that φ is in equilibrium beforeBBN as long as g αβ (cid:38) . × − (MeV /m φ ) . How-ever, as in many other models, this framework is put un-der pressure by the primordial element abundances fromBig Bang Nucleosynthesis (BBN) [11, 32–34]. Incorpo-rating the latest observations, BBN sets a strong con-straint on the eﬀective number of neutrino species  N eﬀ = 2 . ± . . (2)This can be translated into a 2 σ upper bound ∆ N eﬀ < .
42, which severely limits the presence of extra radiationto solve the Hubble problem.In this work, we explore a novel possibility that allowsa large ∆ N eﬀ surpassing the standard BBN constraint in This type of coupling may be connected to the neutrino massgeneration via the Majoron model [25–31], where both scalarand pseudoscalar couplings exist after the spontaneous breakingof lepton number. For singlet Majorons it can be generated bymixing with heavy right-handed Majorana neutrinos in a gaugeinvariant UV completion. a r X i v : . [ h e p - ph ] F e b Eq. (2). In our Majoron-like model given in Eq. (1), with m φ (cid:38) m φ (cid:38) φ ↔ ν + ν reactions after the neutrinos have decoupled from theelectromagnetic plasma at T dec ν ∼ N eﬀ (in-creasing the expansion rate), such that the ﬁnal neutron-to-proton ratio n/p remains almost the same as in thestandard case.After a realistic BBN simulation is performed usingEq. (1), we depict the chi-square function χ as afunction of the scalar mass m φ in the upper panel ofFig. 1 (blue curve). This χ includes the latest mea-surements of the helium-4 mass fraction ( Y P )  andthe deuterium abundance (D / H) , as well as variousnuclear uncertainties. The dotted red curve represents χ obtained simply with Eq. (2), i.e. without anyscalar φ or rise in neutrino temperature, for the given∆ N eﬀ . Parameters with χ > σ level. It can be observed that a ∆ N eﬀ value as largeas 0 . m φ = 1 . χ (cid:39)
2, in contrast to χ (cid:39) N eﬀ value in Eq. (2). In the lower panel ofFig. 1, we also show the preferred baryon-to-photon ratio η ≡ η b × for each m φ . Interestingly, the 1 σ bandaround m φ = 2 MeV matches very well with the indepen-dent determination of η from the CMB ﬁt within themoderately self-interacting neutrino case . In contrast,the standard case, which can be roughly represented by m φ = 10 MeV (cid:29) T dec ν , disagrees with the MI ν value of η by nearly 2 σ .In the following, we will illustrate the idea and resultsin more details. Large extra radiation for BBN . —The improve-ments in the measurement of primordial element abun-dances and cross sections of nuclear reactions have madeBBN an accurate test for physics beyond the StandardModel [35, 41]. The presence of extra radiation duringthe BBN era will accelerate the freeze-out of neutron-proton conversion, resulting in a larger helium-4 abun-dance than the prediction of standard theory. In addi-tion, the abundance of deuterium is extremely sensitiveto the baryon-to-proton ratio η b , leaving BBN essentiallyparameter-free.The model-independent bounds as in Eq. (2) are usu-ally applicable to a “decoupled” ∆ N eﬀ , which is assumedto evolve separately from the Standard Model plasma.For the decoupled ∆ N eﬀ , the main eﬀect is to changethe expansion rate of the Universe, while leaving otheringredients untouched. However, this is not the case if φ tightly couples to neutrinos, such that entropy exchange E x c l ud e dby BB N χ BB N Favored by H measurement Y P + D / H N eff = ± Δ N eff = m ϕ ( MeV ) η MI ν cosmology ( σ ) BB N p r e d i c ti on ( σ ) FIG. 1. The statistical signiﬁcance of BBN χ (upperpanel) and baryon-to-photon ratio η (lower panel) as func-tions of the scalar mass m φ , assuming φ is tightly coupledto all three active neutrinos throughout BBN. In the upperpanel, the solid blue (or dotted red) curve shows χ by fullysimulating nucleosynthesis with the AlterBBN code [39, 40] (oradopting the usual bound N eﬀ = 2 . ± .
27 ). The dashedvertical lines stand for values of ∆ N eﬀ for corresponding m φ .A value ∆ N eﬀ (cid:39) . m φ (cid:38) χ (cid:39)
2, in contrast to the usual BBN bound∆ N eﬀ (cid:46) .
42 at 2 σ level . In the lower panel, the yel-low region signiﬁes the 1 σ allowed range of η predicted byhelium-4 and deuterium abundances for diﬀerent m φ . Theindependent preferred range given by CMB ﬁt  with mod-erately self-interacting neutrinos is shown in the horizontalblue band. It can be noticed that m φ (cid:39) ν . can take place between them. In our case, the argumentbased on ∆ N eﬀ should be taken with caution, and weneed to solve the primordial abundances.Two steps are necessary to derive the light elementabundances. First, the background evolution of vari-ous species ( e ± , γ , ν and φ ) needs to be calculated.Second, we integrate this into a BBN code to numer-ically simulate the synthesis of elements. To calculatethe evolution of background species, we solve the Boltz-mann equations including the weak interactions betweenneutrinos and electrons, so the non-instantaneous decou-pling of neutrinos is taken into account. More detailscan be found in the appendix. We assume that all threegenerations of active neutrinos are in thermal equilib-rium with φ before and during BBN, which holds for g αβ (cid:38) . × − (MeV /m φ ), such that one temperature T ν is adequate to capture the statistical property of the T ν / T γ Entropy transferfrom ϕ m ϕ = Δ N eff = ( decoupled ) Standard n - pfreeze out
10 1 0.1 0.0188.8.131.52.8 T γ ( MeV ) N e ff Ruled out by BBN fora decoupled Δ N eff FIG. 2. The temperature ratio for neutrinos and photons T ν /T γ (upper panel) and N eﬀ (lower panel) with respect tothe photon plasma temperature. For all panels, the bluecurves stand for the case of m φ = 2 . φ isin thermal equilibrium with neutrinos, while the red curvesstand for the case with the decoupled ∆ N eﬀ . The gray curvessigniﬁes the standard case with neutrino-electron decouplingtaken into account. neutrino- φ plasma. This greatly boosts our computationwithout solving discretized distribution functions.In Fig. 2 we show the evolution of temperature ra-tio of neutrino to photon T ν /T γ (upper panel) and N eﬀ (lower panel) as functions of the photon temperature T γ . Two beyond-standard-model scenarios are given: onewith the tightly coupled Majoron-like scalar with mass m φ = 2 . N eﬀ = 0 .
57 (red curves). The standard case with onlythree active neutrinos is shown as gray curves. Note fromthe lower panel that both scenarios are excluded if wesimply adopt the ∆ N eﬀ bound, i.e. if we disregard the ef-fect of φ -interactions on BBN. In the upper panel, for thecase of m φ = 2 . T γ < m φ , the neu-trino plasma receives entropy from the massive φ , and itstemperature is increased by 4 .
6% compared to the stan-dard value. In contrast, for the decoupled ∆ N eﬀ scenariothe ratio T ν /T γ is barely altered. Hence, diﬀerent fromthe decoupled ∆ N eﬀ scenario, there are two eﬀects for thecase m φ = 2 . N eﬀ and higherneutrino temperature T ν . If we assume the entropy from φ is completely transferred to neutrinos, the increasedtemperature can be calculated by using entropy conser-vation. Namely, T ν = ( g ∗ s /g ∗ s ) / T ν = (1 + 6 . T ν ,with g ∗ s ≡ / g ∗ s = 21 / φ decays, respectively.In the realistic case, owing to the weak interactions be-tween neutrinos and electrons, a small part of the entropygoes into the electron-photon plasma. We now investigate these eﬀects on the neutron-to-proton ratio n/p , which is the most important BBN quan-tity before the deuterium bottleneck at T γ (cid:39) .
078 MeVis broken through. For the neutron-proton conversionprocesses where neutrinos appear in the ﬁnal state, e.g.p + e − → n + ν e , the neutrino temperature T ν is rele-vant only through the Pauli blocking factor 1 − f ( p ν e ),which is insensitive to the small change of T ν . Here, f ( p ν e ) stands for the Fermi-Dirac distribution function f ( p ν e ) = 1 / (1 + e p νe /T ν ), where p ν e is the momentum of ν e in the plasma. Thus, we should be concerned aboutonly two processes: n + ν e → p + e − and p + ν e → n + e + .The rates are Γ n ν e = 1 τ n λ m e (cid:90) ∞ d p ν e (cid:113) ( p ν e + Q ) − m e p ν e + Q − ( p νe + Q ) /T γ · p ν e p νe /T ν , (3)Γ p ν e = 1 τ n λ m e (cid:90) ∞ Q + m e d p ν e (cid:113) ( p ν e − Q ) − m e p ν e − Q − ( p νe − Q ) /T γ · p ν e p νe /T ν , (4)where m e is the electron mass, Q ≡ m n − m p (cid:39) .
293 MeV, λ (cid:39) . τ n the neutron lifetime.When T ν < Q , the rate for p + ν e → n + e + is sup-pressed by a Boltzmann factor e − Q/T ν , i.e., only neutri-nos with enough initial energy are kinematically allowedfor the process. In contrast, the neutron-burning processn + ν e → p + e − can take place without energy thresh-old. Hence, after the decoupling of weak interactions at T ν ∼ m e /Q (cid:39) .
15 as a small quantity, the rate inEq. (3) can be well approximated byΓ n ν e (cid:39) τ n λ m e (cid:20) ζ (5)2 T ν + 7 π QT ν + 34 (cid:0) Q − m e (cid:1) ζ (3) T ν (cid:21) . (5)At low neutrino temperatures, the last term will domi-nate, i.e., Γ n ν e ∝ T ν . Consequently, under the small per-turbation of the neutrino temperature δT ν , the rate willbe shifted by δ Γ / Γ n ν e = 3 δT ν /T ν . For the case of m φ =2 . δT ν /T ν = 4 . δ Γ / Γ n ν e (cid:39) . . (cid:46) T ν (cid:46) totn is mainly composed of twoprocesses with similar rates, namely n + ν e → p + e − andn + e − → p + ν e , so we further have δ Γ / Γ totn (cid:39) . m φ = 2 . φ will increase the total conversion rate fromneutrons to protons by almost 6 . N eﬀ . To see that, let us es-timate more precisely the impact of ∆ N eﬀ . Around theBBN era, the Hubble expansion rate is governed by H (cid:39) . √ g ∗ T γ M Pl , (6)where g ∗ = 5 . / · N eﬀ stands for the relativistic de-grees of freedom before the annihilation of electrons, and M Pl = 1 . × GeV for the Planck mass. Note thatthe time scales as t ∝ H − . Hence, under a perturbationof ∆ N eﬀ , the amount of time over a certain temperaturewindow (e.g. from T γ = 1 MeV to T γ = 0 .
078 MeV) willbe changed by δt/t (cid:39) − / · ∆ N eﬀ /g std ∗ with g std ∗ = 10 . N eﬀ = 3 . m φ = 2 . N eﬀ is about 0 .
6, leading to δt/t (cid:39) − . δt/t is negative. Theevolution of neutron number before T γ = 0 . n d t = − Γ totn n + Γ totp p , (7)where n and p are the neutron and proton densities, re-spectively. As has been mentioned before, the conversionrate from proton to neutron Γ totp (cid:39) Γ totn e − Q/T ν is highlysuppressed for T ν < Q (cid:39) .
293 MeV. After the decou-pling of weak interactions at T ν ∼ n/p . Thus, the decreased neutron density in a smallunit time window t is given by Γ totn t n , which is sensitiveto both the perturbations of conversion rate and time(through the expansion rate). As a result, the largerconversion rate with δ Γ / Γ totn (cid:39) .
9% and the larger Hub-ble expansion rate with δH/H (cid:39) .
9% (or δt/t (cid:39) − . AlterBBN code [39, 40] to calculate the light elementabundances, incorporating the background quantities wesolved before (see the appendix). We give in the upperpanel of Fig. 3 the evolution of n/p with respect to thephoton temperature. For illustration, the lower panel in-dicates the diﬀerence of two non-standard scenarios tothe standard one, e.g., n/p | φ − n/p | std . One can clearlynotice the diﬀerent impacts of tightly-coupled φ and thedecoupled ∆ N eﬀ . For the decoupled ∆ N eﬀ , n/p takes alarger value than the standard one (by ∼ . δY p (cid:39) . Y p = 2 n/p (1 + n/p ),which is signiﬁcant given the error of the Y p measure-ment being about 0 .
004 . In contrast, for the case of ascalar m φ = 2 . n/p diﬀers from the standard caseonly by 0 . δY p (cid:39) . n/p increasesinitially due to the larger expansion rate similar to thedecoupled ∆ N eﬀ . Around T γ (cid:39) . n / p m ϕ = Δ N eff = ( decoupled ) Standard T γ ( MeV ) FIG. 3. The neutron-to-proton ratio shown as a function ofthe photon plasma temperature T γ (upper panel). The lowerpanel gives the diﬀerences between the non-standard scenariosand the standard one. The convention is the same as forFig. 2. neutrino temperature starts to accelerate the burning ofneutron, dragging n/p back to the standard value. As aconsequence, Y p is barely altered. This behavior agreesvery well with the previous analytical observations. Preferred parameter space . —Now we explore thepreferred parameter space of the model by setting m φ and g ττ as free parameters, using BBN and other cosmologi-cal observations. The primordial values of light elementabundances can be inferred from the observation of youngastrophysical systems. The mass fraction of helium-4 hasbeen measured to be Y p = 0 . ± . α forests above certain red-shifts. The new recommended deuterium-to-hydrogenabundance ratio reads D / H = (2 . ± . × − .On the other hand, the predicted value of Y p from BBNis dominated by the neutron-to-proton ratio n/p . A stan-dard freeze-out value n/p (cid:39) / N eﬀ (cid:39) Y p (cid:39) .
25. The synthesis of deuterium is morecomplex, depending on both N eﬀ and η b . The remark-able sensitivity of D / H to η b makes it an excellent baryonmeter, especially with the recent update of deuterium-related nuclear rates . In the following, Y p and D / Hwill be used to constrain our model.The mass of the scalar φ cannot be arbitrary in ourscenario. If m φ is too large, the entropy will be mostlyreleased before the neutrino decoupling epoch, and theresulting ∆ N eﬀ is inadequate to explain the Hubble ten-sion. On the other hand, if m φ is too small, there is notenough entropy transfer during the BBN era, and theBBN constraint on ∆ N eﬀ cannot be evaded. This willconﬁne the working range of m φ , as seen in Fig. 1.To fully explore the parameter space of scalar mass m φ and coupling constant g ττ , we incorporate our re-sults into a ﬁt of CMB and large scale structure withself-interacting neutrinos. We note that the observationaldata of CMB and structure formation were initially usedto derive bounds on the secret neutrino interactions [44–53], but later a degeneracy was noticed between the ef-fective neutrino coupling G eﬀ and other cosmological pa-rameters , which can help to resolve the Hubble is-sue. With the Planck 2015 data, the Hubble tension canbe ﬁrmly addressed by a large ∆ N eﬀ along with self-interacting neutrinos . However, the ﬁts based on thelatest Planck 2018 data (speciﬁcally with the high- l po-larization data) show no clear preference for strongly in-teracting neutrinos [7, 8]. The results omitting the high- l polarization data however remain similar to the analysiswith Planck 2015 data. In either case, including the lo-cal measurement of H will always induce a preferencefor large ∆ N eﬀ and non-vanishing G eﬀ , but the overallﬁt with the high- l polarization data of Planck 2018 ispoor.In order to be deﬁnite, we will adopt the results whereonly ν τ moderately couples to φ during the recombina-tion epoch. These include :log ( G eﬀ · MeV ) = − . +1 . − . , N eﬀ = 3 . +0 . − . ,η = 6 . ± . , (8)with the normal neutrino mass ordering. The photon-to-baryon ratio η is converted from Ω b h by using η = 274 Ω b h . For each parameter choice of m φ and g ττ , the total χ is constructed as a combination ofthe BBN one χ and those ﬁtted with the central valueand symmetrized 1 σ error in Eq. (8). The preferred re-gion of parameter space is given in Fig. 4. The dark redregion signiﬁes the preferred parameter space for g ττ at90% conﬁdence level (CL), inside which the dashed curvestands for the 68% CL contour and the star representsthe best-ﬁt point, i.e., m φ = 2 . g ττ = 0 . ν e and ν µ with the scalar. When it comes to the re-combination epoch, to achieve moderate self-interactions G eﬀ ∼ − MeV − , we must have sizable g ττ . Thebound on g ττ from Z decays is shown as the gray bandon the top . On the other hand, to have a higherneutron burning rate, ν e must stay in equilibrium with φ during the BBN era as we assumed in the previousdiscussion, which will impose a lower limit on the cou-pling constant g ee (cid:38) . × − (MeV /m φ ) . The re-quired coupling strength for g ee is depicted in the yellowregion. The lighter yellow regions are excluded by neutri-noless double-beta decay experiments [15, 65] (top-left), K decays  (top-right) and supernova luminosity con- - - - - - m ϕ ( MeV ) C oup li ng ★ ν e i n t h e r m a l e qu ili b r i u m g ττ g ee T oo m u c hh e li u m - N o t e nough Δ N e ff Z decay IceCube FIG. 4. The scalar coupling versus its mass m φ . The redregion is the 90% preferred parameter space for g ττ - m φ bytaking BBN and ﬁts of CMB and the local value of H intoaccount . Note that the results of the moderately self-interacting neutrino mode have been used. Inside the redregion, the dashed curve surrounds the 1 σ region, and thestar in the middle marks the best-ﬁt point. The bound on g ττ from Z decays is shown in the gray band on the top .The required strength of g ee to keep ν e and φ in equilibriumis shown in the yellow region, while various limits to g ee aregiven as lighter yellow regions [11, 15]. The IceCube limits onthe universal couplings are recast as blue curves , whichshould be weakened for the ﬂavor-speciﬁc coupling g ττ . straint [15, 66–70] (bottom), respectively. The presenceof large g ττ coupling will enhance the invisible decay rateof Z , which can be conveniently measured by thenumber of light neutrino species N ν . Our best-ﬁt case m φ = 2 . g ττ = 0 .
07 predicts N ν = 3 . g ττ may lead to a dipin the spectrum of ultra-high energy (UHE) neutrinosobserved at IceCube by scattering oﬀ the relic neutri-nos [54, 58–63]. If we assume the neutrino mass tobe m ν = 0 . E ν = m φ / (2 m ν ) ≈
78 TeV for our best-ﬁt case m φ = 2 . σ parameter space. However, we need tomention that the actual constraints are subject to theneutrino mass spectrum, the neutrino ﬂavor, the modelof sources as well as initial spectrum of UHE neutrinos.For example, in some model where UHE neutrinos aregenerated from decays of dark matter in Milky Way, theconstraints from diﬀuse spectrum do not apply anymore.The constraints in Ref.  will also alter if a diﬀerentspectrum index or ﬂavor-speciﬁc coupling is taken. Concluding remarks . —In this paper, we have ex-plored the role of BBN for the Majoron-like scalar solu-tion in light of the H tension. Note that this work isbased on scalar interactions of Majorana neutrinos, butsimilar or slightly modiﬁed considerations can also bemade for other theories, e.g., complex scalar and vectorinteractions, and even Dirac neutrinos. By numericallysolving the light element abundances, we ﬁnd that a sim-ple Majoron-like scalar with mass (cid:38) MeV can providemoderately self-interaction as well as large ∆ N eﬀ duringthe recombination epoch to relieve the Hubble tension.The widely concerned BBN constraint on ∆ N eﬀ does notapply because it ignores the entropy transfer of a MeV-scale φ that heats up the neutrino bath. The extra radi-ation and the rise in the neutrino temperature are foundto compensate each other, such that the large ∆ N eﬀ iswarranted throughout the BBN era, which should be veryhelpful to address the H tension. Acknowledgement . — GYH would like to thank Kun-Feng Lyu for inspiring discussions. This work is sup-ported by the Alexander von Humboldt Foundation.
In this appendix we explain how the evolution of thebackground is obtained in more details.In the assumption that three ﬂavors of active neutrinosare all tightly coupled to φ , only one temperature T ν issuﬃcient to describe the neutrino- φ plasma. The stateof electron-photon plasma is represented by T γ . To makethe eﬀect of expansion of the Universe explicit, it is conve-nient to introduce the following dimensionless quantitiesin the comoving frame: x ≡ ma, q i ≡ p i · a, z i ≡ T i · a, (cid:101) s i = s i · a , (9)where a is the dimensionful scale factor, x the dimension-less scale factor with m being an arbitrary mass scale (weset m = 1 MeV), q i the comoving momentum for species i , z i the comoving temperature, and (cid:101) s i the comoving en-tropy density. If there is only one massless species inthe Universe, q i , z i and (cid:101) s i will be constant during theexpansion of the Universe.Taking account the weak interactions, the comov-ing entropy density transferred from the electron-photonplasma to the neutrino- φ one can be calculated with [71,72] Hx d (cid:101) s νφ d x = (cid:88) α (cid:90) d q (2 π ) C ν α [ f e , f ν ] ln (cid:18) f ν α − f ν α (cid:19) , (10)where the collision terms C ν α [ f e , f ν ] for six neutrino ﬂa-vors ν α (including antineutrinos) have been widely cal-culated and are available in the literature, and see e.g.,Refs. [72–74]. When the collision terms are vanishing (i.e., no heating from the electron-photon plasma), theentropy in the neutrino- φ plasma is conserved. Also notethat for the neutrino decoupling process, which is not athermal equilibrium process, the total entropy (cid:101) s νφ + (cid:101) s eγ is not preserved in general. The entropy density of theneutrino- φ plasma in terms of the neutrino comovingtemperature reads [71, 72] (cid:101) s νφ = − (cid:88) i = ν,φ (cid:90) d q (2 π ) [ f i ln f i ∓ (1 ± f i ) ln(1 ± f i )] , (11)where the upper and lower signs apply to bosons andfermions, respectively. Here, the distribution functionsfor neutrinos and φ are f ν = 1 / (1 + e q ν /z ν ) and f φ =1 / (1 − e √ q ν + x m φ /m /z ν ), and the sum is over all neutrinospecies and the φ boson. We use the following formulato establish the relation between the variation of entropyand that of the neutrino temperature d z ν / d x :d (cid:101) s νφ d x = ∂ (cid:101) s νφ ∂x + ∂ (cid:101) s νφ ∂z ν · d z ν d x , (12)where ∂ (cid:101) s νφ /∂x and ∂ (cid:101) s νφ /∂z ν can be straightforwardlyobtained with Eq. (11). Some observations on Eq. (12)are very helpful. If we assume there is no entropy trans-ferred from other species, i.e., d (cid:101) s νφ / d x = 0, the varia-tion of comoving temperature of neutrinos d z ν / d x willbe proportional to ∂ (cid:101) s νφ /∂x . For (cid:101) s νφ , the only explicitdependence on x is associated with the mass of φ in f φ ;therefore if m φ (cid:28) T ν , ∂ (cid:101) s νφ /∂x will be negligible and T ν simply follows the red-shift as the Universe expands.The function ∂ (cid:101) s νφ /∂x roughly measures the entropy ﬂowfrom φ to neutrinos.On the other hand, the photon temperature can bederived by utilizing energy conservation x d ρ tot / d x = − ρ tot + P tot ), namely x d T γ d x = − ρ tot + P tot ) − x d T ν d x (cid:16) d ρ γ d T γ + d ρ e d T γ (cid:17) d ρ γ d T γ + d ρ e d T γ , (13)along with m d z i / d x = x d T i / d x + T i for i = ν and γ .The temperature of the photon-electron plasma can alsobe derived by solving the electron collision terms similarto that of neutrinos. But since the energy conservationis a guaranteed result of microscopic processes, they areactually equivalent. By combining Eqs. (10), (12) and(13), we are ready to solve the background quantities forany given initial conditions. ∗ [email protected] † [email protected]  Planck , N. Aghanim et al. , “
Planck 2018 results. VI.Cosmological parameters ,” Astron. Astrophys. (2020) A6, arXiv:1807.06209 . A. G. Riess, S. Casertano, W. Yuan, L. M. Macri, andD. Scolnic, “
Large Magellanic Cloud Cepheid StandardsProvide a 1% Foundation for the Determination of theHubble Constant and Stronger Evidence for Physicsbeyond Λ CDM ,” Astrophys. J. (2019) no. 1, 85, arXiv:1903.07603 . W. L. Freedman et al. , “
The Carnegie-Chicago HubbleProgram. VIII. An Independent Determination of theHubble Constant Based on the Tip of the Red GiantBranch ,” arXiv:1907.05922 . E. Di Valentino et al. , “ Cosmology Intertwined II: TheHubble Constant Tension ,” arXiv:2008.11284 . L. Lancaster, F.-Y. Cyr-Racine, L. Knox, and Z. Pan,“ A tale of two modes: Neutrino free-streaming in theearly universe ,” JCAP (2017) 033, arXiv:1704.06657 . C. D. Kreisch, F.-Y. Cyr-Racine, and O. Dor´e,“ Neutrino puzzle: Anomalies, interactions, andcosmological tensions ,” Phys. Rev. D (2020) no. 12,123505, arXiv:1902.00534 . S. Roy Choudhury, S. Hannestad, and T. Tram,“
Updated constraints on massive neutrinoself-interactions from cosmology in light of the H tension ,” arXiv:2012.07519 . T. Brinckmann, J. H. Chang, and M. LoVerde,“ Self-interacting neutrinos, the Hubble parametertension, and the Cosmic Microwave Background ,” arXiv:2012.11830 . A. Das and S. Ghosh, “ Flavor-speciﬁc InteractionFavours Strong Neutrino Self-coupling ,” arXiv:2011.12315 . A. Mazumdar, S. Mohanty, and P. Parashari, “ Flavourspeciﬁc neutrino self-interaction: H tension andIceCube ,” arXiv:2011.13685 . N. Blinov, K. J. Kelly, G. Z. Krnjaic, and S. D.McDermott, “ Constraining the Self-Interacting NeutrinoInterpretation of the Hubble Tension ,” Phys. Rev. Lett. (2019) no. 19, 191102, arXiv:1905.02727 . K.-F. Lyu, E. Stamou, and L.-T. Wang,“
Self-interacting neutrinos: Solution to Hubble tensionversus experimental constraints ,” Phys. Rev. D (2021) no. 1, 015004, arXiv:2004.10868 . V. Brdar, M. Lindner, S. Vogl, and X.-J. Xu,“
Revisiting neutrino self-interaction constraints from Z and τ decays ,” Phys. Rev. D (2020) no. 11, 115001, arXiv:2003.05339 . F. F. Deppisch, L. Graf, W. Rodejohann, and X.-J. Xu,“ Neutrino Self-Interactions and Double Beta Decay ,”Phys. Rev. D (2020) no. 5, 051701, arXiv:2004.11919 . T. Brune and H. P¨as, “
Massive Majorons andconstraints on the Majoron-neutrino coupling ,” Phys.Rev. D (2019) no. 9, 096005, arXiv:1808.08158 . M. Archidiacono, S. Gariazzo, C. Giunti, S. Hannestad,and T. Tram, “ Sterile neutrino self-interactions: H tension and short-baseline anomalies ,” arXiv:2006.12885 . K. J. Kelly, M. Sen, and Y. Zhang, “ IntimateRelationship Between Sterile Neutrino Dark Matter and ∆ N eﬀ ,” arXiv:2011.02487 . H.-J. He, Y.-Z. Ma, and J. Zheng, “ Resolving Hubble Tension by Self-Interacting Neutrinos with DiracSeesaw ,” JCAP (2020) 003, arXiv:2003.12057 . M. Berbig, S. Jana, and A. Trautner, “ The Hubbletension and a renormalizable model of gauged neutrinoself-interactions ,” Phys. Rev. D (2020) no. 11,115008, arXiv:2004.13039 . O. Seto and Y. Toda, “
Comparing early dark energyand extra radiation solutions to the Hubble tension withBBN ,” arXiv:2101.03740 . F. Arias-Aragon, E. Fernandez-Martinez,M. Gonzalez-Lopez, and L. Merlo, “ Neutrino Massesand Hubble Tension via a Majoron in MFV ,” Eur.Phys. J. C (2021) no. 1, 28, arXiv:2009.01848 . E. Grohs, G. M. Fuller, and M. Sen, “ Consequences ofneutrino self interactions for weak decoupling and bigbang nucleosynthesis ,” JCAP (2020) 001, arXiv:2002.08557 . M. Escudero and S. J. Witte, “ A CMB search for theneutrino mass mechanism and its relation to the Hubbletension ,” Eur. Phys. J. C (2020) no. 4, 294, arXiv:1909.04044 . F. Forastieri, M. Lattanzi, and P. Natoli, “ Cosmologicalconstraints on neutrino self-interactions with a lightmediator ,” Phys. Rev. D (2019) no. 10, 103526, arXiv:1904.07810 . Y. Chikashige, R. N. Mohapatra, and R. D. Peccei,“
Spontaneously Broken Lepton Number andCosmological Constraints on the Neutrino MassSpectrum ,” Phys. Rev. Lett. (1980) 1926. Y. Chikashige, R. N. Mohapatra, and R. D. Peccei, “ AreThere Real Goldstone Bosons Associated with BrokenLepton Number? ,” Phys. Lett. B (1981) 265–268. G. B. Gelmini and M. Roncadelli, “ Left-HandedNeutrino Mass Scale and Spontaneously Broken LeptonNumber ,” Phys. Lett. B (1981) 411–415. K. Choi and A. Santamaria, “ ,” Phys. Lett. B (1991) 504–508. A. Acker, A. Joshipura, and S. Pakvasa, “ A Neutrinodecay model, solar anti-neutrinos and atmosphericneutrinos ,” Phys. Lett. B (1992) 371–375. H. M. Georgi, S. L. Glashow, and S. Nussinov,“
Unconventional Model of Neutrino Masses ,” Nucl.Phys. B (1981) 297–316. J. Schechter and J. W. F. Valle, “
Neutrino Decay andSpontaneous Violation of Lepton Number ,” Phys. Rev.D (1982) 774. G.-y. Huang, T. Ohlsson, and S. Zhou, “ ObservationalConstraints on Secret Neutrino Interactions from BigBang Nucleosynthesis ,” Phys. Rev. D (2018) no. 7,075009, arXiv:1712.04792 . N. Sch¨oneberg, J. Lesgourgues, and D. C. Hooper, “ TheBAO+BBN take on the Hubble tension ,” JCAP (2019) 029, arXiv:1907.11594 . J. Venzor, A. P´erez-Lorenzana, and J. De-Santiago,“ Bounds on neutrino-scalar non-standard interactionsfrom big bang nucleosynthesis ,” arXiv:2009.08104 . C. Pitrou, A. Coc, J.-P. Uzan, and E. Vangioni,“ Precision big bang nucleosynthesis with improvedHelium-4 predictions ,” Phys. Rept. (2018) 1–66, arXiv:1801.08023 .
Neutrino Non-Standard Interactions: A Status Report ,vol. 2. 2019. arXiv:1907.00991 . E. Aver, K. A. Olive, and E. D. Skillman, “
The eﬀects of He I λ ,”JCAP (2015) 011, arXiv:1503.08146 . R. J. Cooke, M. Pettini, and C. C. Steidel, “ OnePercent Determination of the Primordial DeuteriumAbundance ,” Astrophys. J. (2018) no. 2, 102, arXiv:1710.11129 . A. Arbey, “
AlterBBN: A program for calculating theBBN abundances of the elements in alternativecosmologies ,” Comput. Phys. Commun. (2012)1822–1831, arXiv:1106.1363 . A. Arbey, J. Auﬃnger, K. P. Hickerson, and E. S.Jenssen, “
AlterBBN v2: A public code for calculatingBig-Bang nucleosynthesis constraints in alternativecosmologies ,” Comput. Phys. Commun. (2020)106982, arXiv:1806.11095 .
Particle Data Group , P. A. Zyla et al. , “
Review ofParticle Physics ,” PTEP (2020) no. 8, 083C01. S. Weinberg,
Cosmology . 9, 2008. A. Coc, P. Petitjean, J.-P. Uzan, E. Vangioni,P. Descouvemont, C. Iliadis, and R. Longland, “
Newreaction rates for improved primordial D/H calculationand the cosmic evolution of deuterium ,” Phys. Rev. D (2015) no. 12, 123526, arXiv:1511.03843 . S. Hannestad, “ Structure formation with stronglyinteracting neutrinos - Implications for the cosmologicalneutrino mass bound ,” JCAP (2005) 011, arXiv:astro-ph/0411475 . S. Hannestad and G. Raﬀelt, “ Constraining invisibleneutrino decays with the cosmic microwave background ,”Phys. Rev. D (2005) 103514, arXiv:hep-ph/0509278 . N. F. Bell, E. Pierpaoli, and K. Sigurdson,“ Cosmological signatures of interacting neutrinos ,”Phys. Rev. D (2006) 063523, arXiv:astro-ph/0511410 . A. Basboll, O. E. Bjaelde, S. Hannestad, and G. G.Raﬀelt, “ Are cosmological neutrinos free-streaming? ,”Phys. Rev. D (2009) 043512, arXiv:0806.1735 . M. Archidiacono and S. Hannestad, “ Updatedconstraints on non-standard neutrino interactions fromPlanck ,” JCAP (2014) 046, arXiv:1311.3873 . F. Forastieri, M. Lattanzi, and P. Natoli, “ Constraintson secret neutrino interactions after Planck ,” JCAP (2015) 014, arXiv:1504.04999 . F.-Y. Cyr-Racine and K. Sigurdson, “ Limits onNeutrino-Neutrino Scattering in the Early Universe ,”Phys. Rev. D (2014) no. 12, 123533, arXiv:1306.1536 . I. M. Oldengott, C. Rampf, and Y. Y. Y. Wong,“ Boltzmann hierarchy for interacting neutrinos I:formalism ,” JCAP (2015) 016, arXiv:1409.1577 . F. Forastieri, M. Lattanzi, G. Mangano, A. Mirizzi,P. Natoli, and N. Saviano, “ Cosmic microwavebackground constraints on secret interactions amongsterile neutrinos ,” JCAP (2017) 038, arXiv:1704.00626 . I. M. Oldengott, T. Tram, C. Rampf, and Y. Y. Y.Wong, “ Interacting neutrinos in cosmology: exactdescription and constraints ,” JCAP (2017) 027, arXiv:1706.02123 . M. Bustamante, C. Rosenstrøm, S. Shalgar, andI. Tamborra, “ Bounds on secret neutrino interactionsfrom high-energy astrophysical neutrinos ,” Phys. Rev. D (2020) no. 12, 123024, arXiv:2001.04994 .  G. Steigman, “
Primordial Nucleosynthesis in thePrecision Cosmology Era ,” Ann. Rev. Nucl. Part. Sci. (2007) 463–491, arXiv:0712.1100 . Y. Farzan, M. Lindner, W. Rodejohann, and X.-J. Xu,“ Probing neutrino coupling to a light scalar withcoherent neutrino scattering ,” JHEP (2018) 066, arXiv:1802.05171 . G. Arcadi, J. Heeck, F. Heizmann, S. Mertens, F. S.Queiroz, W. Rodejohann, M. Slez´ak, and K. Valerius,“ Tritium beta decay with additional emission of newlight bosons ,” JHEP (2019) 206, arXiv:1811.03530 . K. C. Y. Ng and J. F. Beacom, “ Cosmic neutrinocascades from secret neutrino interactions ,” Phys. Rev.D (2014) no. 6, 065035, arXiv:1404.2288 . [Erratum:Phys.Rev.D 90, 089904 (2014)]. K. Ioka and K. Murase, “ IceCube PeV–EeV neutrinosand secret interactions of neutrinos ,” PTEP (2014) no. 6, 061E01, arXiv:1404.2279 . M. Ibe and K. Kaneta, “
Cosmic neutrino backgroundabsorption line in the neutrino spectrum at IceCube ,”Phys. Rev. D (2014) no. 5, 053011, arXiv:1407.2848 . A. Kamada and H.-B. Yu, “ Coherent Propagation ofPeV Neutrinos and the Dip in the Neutrino Spectrum atIceCube ,” Phys. Rev. D (2015) no. 11, 113004, arXiv:1504.00711 . I. M. Shoemaker and K. Murase, “ Probing BSMNeutrino Physics with Flavor and Spectral Distortions:Prospects for Future High-Energy Neutrino Telescopes ,”Phys. Rev. D (2016) no. 8, 085004, arXiv:1512.07228 . A. DiFranzo and D. Hooper, “ Searching for MeV-ScaleGauge Bosons with IceCube ,” Phys. Rev. D (2015)no. 9, 095007, arXiv:1507.03015 . S. Shalgar, I. Tamborra, and M. Bustamante,“ Core-collapse supernovae stymie secret neutrinointeractions ,” arXiv:1912.09115 . K. Blum, Y. Nir, and M. Shavit, “ Neutrinolessdouble-beta decay with massive scalar emission ,” Phys.Lett. B (2018) 354–361, arXiv:1802.08019 . E. W. Kolb and M. S. Turner, “
Supernova SN 1987aand the Secret Interactions of Neutrinos ,” Phys. Rev. D (1987) 2895. R. V. Konoplich and M. Y. Khlopov, “ Constraints ontriplet Majoron model due to observations of neutrinosfrom stellar collapse ,” Sov. J. Nucl. Phys. (1988)565–566. Y. Farzan, “ Bounds on the coupling of the Majoron tolight neutrinos from supernova cooling ,” Phys. Rev. D (2003) 073015, arXiv:hep-ph/0211375 . S. Zhou, “ Comment on astrophysical consequences of aneutrinophilic 2HDM ,” Phys. Rev. D (2011) 038701, arXiv:1106.3880 . L. Heurtier and Y. Zhang, “ Supernova Constraints onMassive (Pseudo)Scalar Coupling to Neutrinos ,” JCAP (2017) 042, arXiv:1609.05882 . J. Bernstein, KINETIC THEORY IN THEEXPANDING UNIVERSE . Cambridge Monographs onMathematical Physics. Cambridge University Press,Cambridge, U.K., 1988. E. Grohs, G. M. Fuller, C. T. Kishimoto, M. W. Paris,and A. Vlasenko, “
Neutrino energy transport in weakdecoupling and big bang nucleosynthesis ,” Phys. Rev. D (2016) no. 8, 083522, arXiv:1512.02205 .  S. Hannestad and J. Madsen, “ Neutrino decoupling inthe early universe ,” Phys. Rev. D (1995) 1764–1769, arXiv:astro-ph/9506015 . A. D. Dolgov, S. H. Hansen, and D. V. Semikoz, “ Nonequilibrium corrections to the spectra of masslessneutrinos in the early universe ,” Nucl. Phys. B (1997) 426–444, arXiv:hep-ph/9703315arXiv:hep-ph/9703315