Solving the Hubble tension without spoiling Big Bang Nucleosynthesis

The Hubble parameter inferred from cosmic microwave background observations is consistently lower than that from local measurements, which could hint towards new physics. Solutions to the Hubble tension typically require a sizable amount of extra radiation $\Delta N^{}_{\rm eff}$ during recombination. However, the amount of $\Delta N^{}_{\rm eff}$ 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 problem by introducing neutrino self-interactions via a simple Majoron-like coupling. The scalar is slightly heavier than $1~{\rm MeV}$ and allowed to be fully thermalized throughout the BBN era. The rise of neutrino temperature due to the entropy transfer via $\phi \to \nu\overline{\nu}$ reactions compensates the effect of a large $\Delta N^{}_{\rm eff}$ on BBN. Values of $\Delta N^{}_{\rm eff}$ as large as $0.7$ are in this case compatible with BBN. We perform a fit to the parameter space of the model.

The Hubble parameter inferred from cosmic microwave background observations is consistently lower than that from local measurements, which could hint towards new physics. Solutions to the Hubble tension typically require a sizable amount of extra radiation ∆N eff during recombination. However, the amount of ∆N eff 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 problem by introducing neutrino self-interactions via a simple Majoron-like coupling. The scalar is slightly heavier than 1 MeV and allowed to be fully thermalized throughout the BBN era. The rise of neutrino temperature due to the entropy transfer via φ → νν reactions compensates the effect of a large ∆N eff on BBN. Values of ∆N eff as large as 0.7 are in this case compatible with BBN. We perform a fit to the parameter space of the model.
Introduction.-The Hubble parameter inferred from the Planck observations of the cosmic microwave background (CMB), H 0 = 67.4 ± 0.5 km/s/Mpc [1], is in tension with that of local measurements at low red-shifts. To be specific, the result from the Hubble Space Telescope (HST) by observing the Milky Way Cepheids is H 0 = 74.03±1.42 km/s/Mpc [2], which exceeds the value of Planck experiment by a 4.4σ significance. Combining the HST result with a later independent determination [3] yields a lower value of H 0 = 72.26 ± 1.19 km/s/Mpc, but the tension still persists at 3.7σ level.
The tension for the Hubble parameter could suggest the existence of new physics beyond the Standard Model or beyond the ΛCDM framework [4]. There is a strong positive correlation between H 0 and an extra radiation, ∆N eff = N eff −3.046, in the early Universe. Hence, by increasing N eff during recombination one can lift the Hubble parameter. However, increasing N eff delays matterradiation equality and modifies the CMB power spectrum. This, in turn, can be compensated by introducing non-standard neutrino self-interactions during recombination [5,6]. Thus, a successful particle physics model to explain the Hubble tension needs to provide a significant amount of ∆N eff in the early Universe, as well as "secret" neutrino interactions. In the original fits with Planck 2015 data, two modes with self-interacting neutrinos of the form G eff νννν are identified, which are given in Table I [6]. SIν (MIν) stands for the strongly (moderately) interacting neutrino mode; η 10 ≡ η b × 10 10 represents the baryon-to-photon ratio. These results are updated with the Planck 2018 data in Refs. [7][8][9][10].
In our model, the scalar particle φ with a mass m φ increases N eff . As long as the coupling in Eq. (1) is strong enough, φ will be in thermal equilibrium with the neutrino plasma, contributing to extra radiation by ∆N eff = 1/2 · 8/7 0.57 for m φ T ν , where T ν is the plasma temperature. Note that φ is in equilibrium before BBN as long as g αβ 2.2 × 10 −10 (MeV/m φ ) [32]. However, as in many other models, this framework is put under pressure by the primordial element abundances from Big Bang Nucleosynthesis (BBN) [11,[32][33][34]. Incorporating the latest observations, BBN sets a strong constraint on the effective number of neutrino species [35] This can be translated into a 2σ upper bound ∆N eff < 0.42, which severely limits the presence of extra radiation to solve the Hubble problem.
In this work, we explore a novel possibility that allows a large ∆N eff surpassing the standard BBN constraint in Eq. (2). In our Majoron-like model given in Eq. (1), with m φ 1 MeV, the scalar particle can stay safely in thermal equilibrium throughout the epoch of BBN, in contrast to concerns in the literature [11,33,36]. Namely, since m φ 1 MeV, the neutrino temperature will increase with respect to the photon one due to φ ↔ ν + ν reactions after the neutrinos have decoupled from the electromagnetic plasma at T dec ν ∼ 1 MeV. The rise in the neutrino temperature (increasing the neutron burning rate) will cancel the effect caused by a larger N eff (increasing the expansion rate), such that the final neutronto-proton ratio n/p remains almost the same as in the standard case.
After a realistic BBN simulation is performed using Eq. (1), we depict the chi-square function χ 2 BBN as a function of the scalar mass m φ in the upper panel of Fig. 1 (blue curve). This χ 2 BBN includes the latest measurements of the helium-4 mass fraction (Y P ) [37] and the deuterium abundance (D/H) [38], as well as various nuclear uncertainties. The dotted red curve represents χ 2 BBN obtained simply with Eq. (2), i.e. without any scalar φ or rise in neutrino temperature, for the given ∆N eff . Parameters with χ 2 BBN > 4 are ruled out at 2σ level. It can be observed that a ∆N eff value as large as 0.7 is allowed for m φ = 1.8 MeV without spoiling BBN, i.e., χ 2 BBN 2, in contrast to χ 2 BBN 9 using simply the N eff value in Eq. (2). In the lower panel of Fig. 1, we also show the preferred baryon-to-photon ratio η 10 ≡ η b × 10 10 for each m φ . Interestingly, the 1σ band around m φ = 2 MeV matches very well with the independent determination of η 10 from the CMB fit within the moderately self-interacting neutrino case [6]. In contrast, the standard case, which can be roughly represented by m φ = 10 MeV T dec ν , disagrees with the MIν value of η 10 by nearly 2σ.
In the following, we will illustrate the idea and results in more details.
Large extra radiation for BBN.-The improvements in the measurement of primordial element abundances and cross sections of nuclear reactions have made BBN an accurate test for physics beyond the Standard Model [35,41]. The presence of extra radiation during the BBN era will accelerate the freeze-out of neutronproton conversion, resulting in a larger helium-4 abundance than the prediction of standard theory. In addition, the abundance of deuterium is extremely sensitive to the baryon-to-proton ratio η b , leaving BBN essentially parameter-free.
The model-independent bounds as in Eq. (2) are usually applicable to a "decoupled" ∆N eff , which is assumed to evolve separately from the Standard Model plasma. For the decoupled ∆N eff , the main effect is to change the expansion rate of the Universe, while leaving other ingredients untouched. However, this is not the case if φ tightly couples to neutrinos, such that entropy exchange  [35]. In the lower panel, the yellow region signifies the 1σ allowed range of η 10 predicted by helium-4 and deuterium abundances for different m φ . The independent preferred range given by CMB fit [6] with moderately self-interacting neutrinos is shown in the horizontal blue band. It can be noticed that m φ 2 MeV provides excellent fits to both BBN and CMB with MIν.
can take place between them. In our case, the argument based on ∆N eff should be taken with caution, and we need to solve the primordial abundances.
Two steps are necessary to derive the light element abundances. First, the background evolution of various species (e ± , γ, ν and φ) needs to be calculated. Second, we integrate this into a BBN code to numerically simulate the synthesis of elements. To calculate the evolution of background species, we solve the Boltzmann equations including the weak interactions between neutrinos and electrons, so the non-instantaneous decoupling of neutrinos is taken into account. More details can be found in the appendix. We assume that all three generations of active neutrinos are in thermal equilibrium with φ before and during BBN, which holds for g αβ 2.2 × 10 −10 (MeV/m φ ), such that one temperature T ν is adequate to capture the statistical property of the neutrino-φ plasma. This greatly boosts our computation without solving discretized distribution functions.
In Fig. 2 we show the evolution of temperature ratio of neutrino to photon T ν /T γ (upper panel) and N eff (lower panel) as functions of the photon temperature T γ . Two beyond-standard-model scenarios are given: one with the tightly coupled Majoron-like scalar with mass m φ = 2.4 MeV (blue curves), and one with the decoupled ∆N eff = 0.57 (red curves). The standard case with only three active neutrinos is shown as gray curves. Note from the lower panel that both scenarios are excluded if we simply adopt the ∆N eff bound, i.e. if we disregard the effect of φ-interactions on BBN. In the upper panel, for the case of m φ = 2.4 MeV, shortly after T γ < m φ , the neutrino plasma receives entropy from the massive φ, and its temperature is increased by 4.6% compared to the standard value. In contrast, for the decoupled ∆N eff scenario the ratio T ν /T γ is barely altered. Hence, different from the decoupled ∆N eff scenario, there are two effects for the case m φ = 2.4 MeV: extra radiation ∆N eff and higher neutrino temperature T ν . If we assume the entropy from φ is completely transferred to neutrinos, the increased temperature can be calculated by using entropy conservation. Namely, T 0 ν = (g * s /g 0 * s ) 1/3 T ν = (1 + 6.0%)T ν , with g * s ≡ 25/4 and g 0 * s = 21/4 being the entropy degrees of freedom before and after φ decays, respectively. In the realistic case, owing to the weak interactions between neutrinos and electrons, a small part of the entropy goes into the electron-photon plasma.
We now investigate these effects on the neutron-toproton ratio n/p, which is the most important BBN quantity before the deuterium bottleneck at T γ 0.078 MeV is broken through. For the neutron-proton conversion processes where neutrinos appear in the final state, e.g. p + e − → n + ν e , the neutrino temperature T ν is relevant 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 about only two processes: n + ν e → p + e − and p + ν e → n + e + . The rates are [42] Γ nν e = 1 τ n λ m 5 where m e is the electron mass, Q ≡ m n − m p 1.293 MeV, λ 1.636, and τ n the neutron lifetime. When T ν < Q, the rate for p + ν e → n + e + is suppressed by a Boltzmann factor e −Q/T ν , i.e., only neutrinos with enough initial energy are kinematically allowed for the process. In contrast, the neutron-burning process n + ν e → p + e − can take place without energy threshold. Hence, after the decoupling of weak interactions at T ν ∼ 1 MeV, a higher neutrino temperature compared to the standard case will result in a larger neutron burning rate. By ignoring the electron distribution function in the Pauli blocking factor and expanding the square-root by taking m e /Q 0.15 as a small quantity, the rate in Eq. (3) can be well approximated by At low neutrino temperatures, the last term will dominate, i.e., Γ nν e ∝ T 3 ν . Consequently, under the small perturbation of the neutrino temperature δT ν , the rate will be shifted by δΓ/Γ nν e = 3 δT ν /T ν . For the case of m φ = 2.4 MeV in Fig. 2, the relative temperature shift is about δT ν /T ν = 4.6%, so we have δΓ/Γ nν e 13.8%. During the temperature window 0.2 MeV T ν 1 MeV, the total neutron conversion rate Γ tot n is mainly composed of two processes with similar rates, namely n + ν e → p + e − and n + e − → p + ν e , so we further have δΓ/Γ tot n 6.9%. To conclude, for the case with m φ = 2.4 MeV, the change of neutrino temperature induced by the entropy transfer from φ will increase the total conversion rate from neutrons to protons by almost 6.9%.
The above result will compensate the larger expansion rate caused by a positive ∆N eff . To see that, let us estimate more precisely the impact of ∆N eff . Around the BBN era, the Hubble expansion rate is governed by where g * = 5.5 + 7/4 · N eff stands for the relativistic degrees of freedom before the annihilation of electrons, and M Pl = 1.221 × 10 19 GeV for the Planck mass. Note that the time scales as t ∝ H −1 . Hence, under a perturbation of ∆N eff , the amount of time over a certain temperature window (e.g. from T γ = 1 MeV to T γ = 0.078 MeV) will be changed by δt/t −7/8·∆N eff /g std * with g std * = 10.75 being the degrees of freedom with N eff = 3.046. For our case of m φ = 2.4 MeV, ∆N eff is about 0.6, leading to δt/t −4.9%. It is crucial that δt/t is negative. The evolution of neutron number before T γ = 0.1 MeV is described by dn dt = −Γ tot n n + Γ tot p p , where n and p are the neutron and proton densities, respectively. As has been mentioned before, the conversion rate from proton to neutron Γ tot p Γ tot n e −Q/T ν is highly suppressed for T ν < Q 1.293 MeV. After the decoupling of weak interactions at T ν ∼ 1 MeV, the conversion from neutrons to protons will dominate the evolution of n/p. Thus, the decreased neutron density in a small unit time window t is given by Γ tot n t n, which is sensitive to both the perturbations of conversion rate and time (through the expansion rate). As a result, the larger conversion rate with δΓ/Γ tot n 6.9% and the larger Hubble expansion rate with δH/H 4.9% (or δt/t −4.9%) will compensate each other.
Having some analytical understanding, we adopt the AlterBBN code [39,40] to calculate the light element abundances, incorporating the background quantities we solved before (see the appendix). We give in the upper panel of Fig. 3 the evolution of n/p with respect to the photon temperature. For illustration, the lower panel indicates the difference of two non-standard scenarios to the standard one, e.g., n/p| φ − n/p| std . One can clearly notice the different impacts of tightly-coupled φ and the decoupled ∆N eff . For the decoupled ∆N eff , n/p takes a larger value than the standard one (by ∼ 0.005) due to the higher expansion rate. This leads to a larger helium-4 abundance by δY p 0.0076, with Y p = 2n/p(1 + n/p), which is significant given the error of the Y p measurement being about 0.004 [37]. In contrast, for the case of a scalar m φ = 2.4 MeV, n/p differs from the standard case only by 0.001, resulting in a negligible change of helium-4 abundance δY p 0.0015. In this scenario, n/p increases initially due to the larger expansion rate similar to the decoupled ∆N eff . Around T γ 0.9 MeV, the increasing neutrino temperature starts to accelerate the burning of neutron, dragging n/p back to the standard value. As a consequence, Y p is barely altered. This behavior agrees very well with the previous analytical observations. Preferred parameter space.-Now we explore the preferred parameter space of the model by setting m φ and g τ τ as free parameters, using BBN and other cosmological observations. The primordial values of light element abundances can be inferred from the observation of young astrophysical systems. The mass fraction of helium-4 has been measured to be Y p = 0.2449 ± 0.0040 by observing the emission spectrum of low-metallicity compact blue galaxies [37]. In addition, the deuterium abundance with a much lower value was derived by observing the absorption spectrum of Lyman-α forests above certain redshifts. The new recommended deuterium-to-hydrogen abundance ratio reads D/H = (2.527±0.030)×10 −5 [38]. On the other hand, the predicted value of Y p from BBN is dominated by the neutron-to-proton ratio n/p. A standard freeze-out value n/p 1/7 with N eff 3 will give rise to Y p 0.25. The synthesis of deuterium is more complex, depending on both N eff and η b . The remarkable sensitivity of D/H to η b makes it an excellent baryon meter, especially with the recent update of deuteriumrelated nuclear rates [43]. In the following, Y p and D/H will be used to constrain our model.
The mass of the scalar φ cannot be arbitrary in our scenario. If m φ is too large, the entropy will be mostly released before the neutrino decoupling epoch, and the resulting ∆N eff is inadequate to explain the Hubble tension. On the other hand, if m φ is too small, there is not enough entropy transfer during the BBN era, and the BBN constraint on ∆N eff cannot be evaded. This will confine 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 results into a fit of CMB and large scale structure with self-interacting neutrinos. We note that the observational data of CMB and structure formation were initially used to derive bounds on the secret neutrino interactions [44][45][46][47][48][49][50][51][52][53], but later a degeneracy was noticed between the effective neutrino coupling G eff and other cosmological parameters [5], which can help to resolve the Hubble issue. With the Planck 2015 data, the Hubble tension can be firmly addressed by a large ∆N eff along with selfinteracting neutrinos [6]. However, the fits based on the latest Planck 2018 data (specifically with the high-l polarization data) show no clear preference for strongly interacting neutrinos [7,8]. The results omitting the high-l polarization data however remain similar to the analysis with Planck 2015 data. In either case, including the local measurement of H 0 will always induce a preference for large ∆N eff and non-vanishing G eff , but the overall fit with the high-l polarization data of Planck 2018 is poor.
In order to be definite, we will adopt the results where only ν τ moderately couples to φ during the recombination epoch. These include [10]: with the normal neutrino mass ordering. The photonto-baryon ratio η 10 is converted from Ω b h 2 by using η 10 = 274 Ω b h 2 [55]. For each parameter choice of m φ and g τ τ , the total χ 2 is constructed as a combination of the BBN one χ 2 BBN and those fitted with the central value and symmetrized 1σ error in Eq. (8). The preferred region of parameter space is given in Fig. 4. The dark red region signifies the preferred parameter space for g τ τ at 90% confidence level (CL), inside which the dashed curve stands for the 68% CL contour and the star represents the best-fit point, i.e., m φ = 2.8 MeV and g τ τ = 0.07.
Laboratory and astrophysical searches set stringent upper limits on the secret neutrino interactions [12-15, 54, 56-64], especially for the coupling strengths of ν e and ν µ with the scalar. When it comes to the recombination epoch, to achieve moderate self-interactions G eff ∼ 10 −3 MeV −2 , we must have sizable g τ τ . The bound on g τ τ from Z decays is shown as the gray band on the top [13]. On the other hand, to have a higher neutron burning rate, ν e must stay in equilibrium with φ during the BBN era as we assumed in the previous discussion, which will impose a lower limit on the coupling constant g ee 2.2 × 10 −10 (MeV/m φ ) [32]. The required coupling strength for g ee is depicted in the yellow region. The lighter yellow regions are excluded by neutrinoless double-beta decay experiments [15,65] (top-left), K decays [11] (top-right) and supernova luminosity con- The red region is the 90% preferred parameter space for g τ τ -m φ by taking BBN and fits of CMB and the local value of H 0 into account [10]. Note that the results of the moderately selfinteracting neutrino mode have been used. Inside the red region, the dashed curve surrounds the 1σ region, and the star in the middle marks the best-fit point. The bound on g τ τ from Z decays is shown in the gray band on the top [13]. The required strength of g ee to keep ν e and φ in equilibrium is shown in the yellow region, while various limits to g ee are given as lighter yellow regions [11,15]. The IceCube limits on the universal couplings are recast as blue curves [54], which should be weakened for the flavor-specific coupling g τ τ .
straint [15,[66][67][68][69][70] (bottom), respectively. The presence of large g τ τ coupling will enhance the invisible decay rate of Z [13], which can be conveniently measured by the number of light neutrino species N ν . Our best-fit case m φ = 2.8 MeV and g τ τ = 0.07 predicts N ν = 3.0012 [13]. The region of our interest for g τ τ may lead to a dip in the spectrum of ultra-high energy (UHE) neutrinos observed at IceCube by scattering off the relic neutrinos [54,[58][59][60][61][62][63]. If we assume the neutrino mass to be m ν = 0.1 eV, the resonant-scattering dip should be around E ν = m 2 φ /(2m ν ) ≈ 78 TeV for our best-fit case m φ = 2.8 MeV. The absence of the dip at IceCube will place a constraint on our preferred parameter space [54], recast as dotted blue curves in Fig. 4, which has covered part of ours 1σ parameter space. However, we need to mention that the actual constraints are subject to the neutrino mass spectrum, the neutrino flavor, the model of sources as well as initial spectrum of UHE neutrinos. For example, in some model where UHE neutrinos are generated from decays of dark matter in Milky Way, the constraints from diffuse spectrum do not apply anymore. The constraints in Ref. [54] will also alter if a different spectrum index or flavor-specific coupling is taken.
Concluding remarks.-In this paper, we have explored the role of BBN for the Majoron-like scalar solution in light of the H 0 tension. Note that this work is based on scalar interactions of Majorana neutrinos, but similar or slightly modified considerations can also be made for other theories, e.g., complex scalar and vector interactions, and even Dirac neutrinos. By numerically solving the light element abundances, we find that a simple Majoron-like scalar with mass MeV can provide moderately self-interaction as well as large ∆N eff during the recombination epoch to relieve the Hubble tension. The widely concerned BBN constraint on ∆N eff does not apply because it ignores the entropy transfer of a MeVscale φ that heats up the neutrino bath. The extra radiation and the rise in the neutrino temperature are found to compensate each other, such that the large ∆N eff is warranted throughout the BBN era, which should be very helpful to address the H 0 tension.
Acknowledgement.-GYH would like to thank Kun-Feng Lyu for inspiring discussions. This work is supported by the Alexander von Humboldt Foundation.

Appendix
In this appendix we explain how the evolution of the background is obtained in more details.
In the assumption that three flavors of active neutrinos are all tightly coupled to φ, only one temperature T ν is sufficient to describe the neutrino-φ plasma. The state of electron-photon plasma is represented by T γ . To make the effect of expansion of the Universe explicit, it is convenient to introduce the following dimensionless quantities in the comoving frame: where a is the dimensionful scale factor, x the dimensionless scale factor with m being an arbitrary mass scale (we set m = 1 MeV), q i the comoving momentum for species i, z i the comoving temperature, and s i the comoving entropy density. If there is only one massless species in the Universe, q i , z i and s i will be constant during the expansion of the Universe.
Taking account the weak interactions, the comoving entropy density transferred from the electron-photon plasma to the neutrino-φ one can be calculated with [71,72] Hx where the collision terms C να [f e , f ν ] for six neutrino flavors ν α (including antineutrinos) have been widely calculated and are available in the literature, and see e.g., Refs. [72][73][74]. When the collision terms are vanishing (i.e., no heating from the electron-photon plasma), the entropy in the neutrino-φ plasma is conserved. Also note that for the neutrino decoupling process, which is not a thermal equilibrium process, the total entropy s νφ + s eγ is not preserved in general. The entropy density of the neutrino-φ plasma in terms of the neutrino comoving temperature reads [71,72] where the upper and lower signs apply to bosons and fermions, respectively. Here, the distribution functions for neutrinos and φ are f ν = 1/(1 + e q ν /z ν ) and f φ = 1/(1−e √ q 2 ν +x 2 m 2 φ /m 2 /z ν ), and the sum is over all neutrino species and the φ boson. We use the following formula to establish the relation between the variation of entropy and that of the neutrino temperature dz ν /dx: where ∂ s νφ /∂x and ∂ s νφ /∂z ν can be straightforwardly obtained with Eq. (11). Some observations on Eq. (12) are very helpful. If we assume there is no entropy transferred from other species, i.e., d s νφ /dx = 0, the variation of comoving temperature of neutrinos dz ν /dx will be proportional to ∂ s νφ /∂x. For s νφ , the only explicit dependence on x is associated with the mass of φ in f φ ; therefore if m φ T ν , ∂ s νφ /∂x will be negligible and T ν simply follows the red-shift as the Universe expands. The function ∂ s νφ /∂x roughly measures the entropy flow from φ to neutrinos.
On the other hand, the photon temperature can be derived by utilizing energy conservation xdρ tot /dx = −3(ρ tot + P tot ), namely along with m dz i /dx = x dT i /dx + T i for i = ν and γ.
The temperature of the photon-electron plasma can also be derived by solving the electron collision terms similar to that of neutrinos. But since the energy conservation is a guaranteed result of microscopic processes, they are actually equivalent. By combining Eqs. (10), (12) and (13), we are ready to solve the background quantities for any given initial conditions.