Precision CMB constraints on eV-scale bosons coupled to neutrinos

The cosmic microwave background (CMB) has proven to be an invaluable tool for studying the properties and interactions of neutrinos, providing insight not only into the sum of neutrino masses but also the free streaming nature of neutrinos prior to recombination. The CMB is a particularly powerful probe of new eV-scale bosons interacting with neutrinos, as these particles can thermalize with neutrinos via the inverse decay process, νν¯→X\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu \bar{\nu } \rightarrow X$$\end{document}, and suppress neutrino free streaming near recombination – even for couplings as small as λν∼O(10-13)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _{\nu } \sim {\mathcal {O}}(10^{-13})$$\end{document}. Here, we revisit CMB constraints on such bosons, improving upon a number of approximations previously adopted in the literature and generalizing the constraints to a broader class of models. This includes scenarios in which the boson is either spin-0 or spin-1, the number of interacting neutrinos is either Nint=1,2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_{\textrm{int}} = 1,2 $$\end{document} or 3, and the case in which a primordial abundance of the species is present. We apply these bounds to well-motivated models, such as the singlet majoron model or a light U(1)Lμ-Lτ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$U(1)_{L_{\mu }-L_{\tau }}$$\end{document} gauge boson, and find that they represent the leading constraints for masses mX∼1eV\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_X\sim 1\, {\textrm{eV}}$$\end{document}. Finally, we revisit the extent to which neutrino-philic bosons can ameliorate the Hubble tension, and find that recent improvements in the understanding of how such bosons damp neutrino free streaming reduces the previously found success of this proposal.


I. INTRODUCTION
Neutrinos always comprise a sizable fraction of the energy density in the Universe.In particular, prior to matter-radiation equality they represent ∼ 40% of the energy budget.Neutrinos are also the only species with a sizable anisotropic stress -a consequence of their decoupling from the thermal plasma at T ∼ 2 MeV.Collectively, these facts imply that neutrino free streaming plays an important role in the evolution of the gravitational potentials responsible for sourcing the CMB anisotropies [1][2][3].Current observations of the CMB by the Planck satellite [4][5][6] are compatible with the standard picture in which neutrinos are free streaming at redshifts 2000 ≲ z ≲ 10 5 [7] (corresponding to temperatures 0.5 eV ≲ T γ ≲ 25 eV), implying these observations can be used to stringently constrain the existence of new light particles coupled to the neutrino sector.
The goal of this work is to perform a precision study of the impact of eV−scale neutrino-philic bosons on the CMB, improving upon previous analyses which relied on numerous simplified approximations [33][34][35], and extending the results of these analyses to the more general class of light neutrino-philic bosons.The primary improvements of this work are three-fold.First, we have incorporated the background thermodynamic evolution of neutrinos and the neutrino-philic bosons in the cosmological Boltzmann code CLASS [46,47].This allows us to solve for the thermodynamics on the fly, with precision and speed which allows a full Bayesian analysis of Planck legacy data 1 .Next, we incorporated a refined computa- 1 Our modified version of CLASS is available on github .The FIG. 1. Parameter space for neutrino interactions with a scalar (left panel ) and vector (right panel ) boson X with mass mX .The bounds are interpreted within the singlet majoron model, where λν = mν /vL and for a light U (1)L µ−Lτ gauge boson, for which λν ≃ gµ−τ , respectively.An analysis of Planck legacy data excludes blue regions with 3σ confidence.Grey regions represent current cosmological, astrophysical and laboratory constraints, see Section V for details.In pink we indicate constraints coming from the out-of-equilibrium decay of the new X boson which apply if a primordial abundance was generated before BBN.We also indicate the region of parameter space which will be tested by the Simons Observatory.In particular, the region above the purple dashed-dotted line will be tested because the thermalization of the X boson leads to an observable excess of ∆N eff ≥ 0.1.Finally, we also highlight in red the best fit region of parameter space for the scenario of the X boson being of scalar type and interacting with one neutrino family, Nint = 1.This region is of particular interest because it indicates that non-trivial neutrino interactions are statistically slightly preferred over ΛCDM.
tion of the collision term [30,31] which damps the neutrino free streaming less efficiently than assumed in previous studies [33][34][35].Finally, we generalize the analysis to arbitrary number of interacting neutrino species, include the possibility of both vector and scalar bosons and the possibility of having a primordial abundance such bosons.
In general, we find that the CMB can robustly constrain the existence of eV−scale neutrino-philic bosons with couplings on the order of λ ν ∼ O(10 −13 ).The value of this coupling roughly corresponds to the new bosonic particles having a lifetime shorter than the age of the Universe at recombination, Γ X ∼ λ 2 ν m X /(8π) ≲ H(z rec ).These bounds play an important role in testing a variety of well-motivated high-energy theories, such as the singlet majoron model (where these observations are testing scales of lepton number breaking as high as ∼ 1 TeV), and the U (1) Lµ−Lτ extension of the Standard Model.The main results of our study are highlighted in Figure 1, which display the 3σ constraint on the coupling of the majoron and U (1) Lµ−Lτ gauge boson, respectively.In the case of the majoron, we also highlight a region of parameter space that is favoured by Planck legacy data at the ∼ 1σ level.
The reminder of this work is structured as follows.First, in Section II we briefly introduce and motivate the particle physics models that we consider.In Section III, we present the formalism behind our work.In particular, we describe how we treat the thermodynamic evolution of equations for the evolution of the temperature and chemical potentials should be easily generalizable to other scenarios involving Beyond the Standard Model (BSM) physics.
the Universe in the presence of eV−scale neutrino-philic bosons, including how the dynamics are implemented at the level of both the background and perturbations.In Section IV we present the constraints we derive on the couplings between neutrinos and eV−scale bosons.We also include a quantitative discussion about the ability of these models to solve or ameliorate the Hubble tension, showing that the new collision term strongly suppresses the previous success of this model identified in [33][34][35].Finally, in Section VI we present a summary of our results and outline our conclusions.For completeness, we provide in the appendices I and II further information on the formalism and details on the modified cosmological history.

II. PARTICLE PHYSICS MODELS
Effective Interactions: We will consider an effective coupling between neutrinos and a light bosonic mediator X and we will study two cases, one where the mediator is a pseudoscalar X = ϕ and one where it is a vector X = Z ′ .We will work after electroweak symmetry breaking and in the active neutrino mass basis.The effective Lagrangians describing these interactions are: where λ ν are dimensionless coupling constants and where a = 1 for Majorana neutrinos and a = √ 2 for Dirac neutrinos.
Given these interactions, the scalar and vector boson partial decay rate into a pair of massive neutrinos are given by: Mapping to concrete models: These effective Lagrangians have a direct interpretation in terms of well motivated BSM scenarios.For example, Eq. ( 1) is the effective interaction generated in the famous singlet majoron model [36] with X = ϕ identified as the majoron and with λ ν = m ν /v L , where v L is the scale at which the global U (1) L symmetry is spontaneously broken.In particular, in this model the coupling between massive neutrinos and the majoron is diagonal up to small corrections [37].The vector interactions in Eq. ( 2) also effectively describe new interactions of neutrinos in many BSM constructions.Typically, in the vector case the interaction arises by the gauging of lepton number family symmetries, and as such, the interaction is non-diagonal in the neutrino mass basis [40,41].However, in such cases all massive neutrinos couple to the X boson, and the couplings in the mass and flavor basis are simply related by a PMNS rotation.As an example, we can consider the case of a light U (1) Lµ−Lτ gauge boson; here, the coupling λ ν is intimately related to the U (1) Lµ−Lτ gauge coupling, λ ν ≃ g µ−τ -see Ref. [48] for the precise mapping.
Scenarios Considered: We will consider several scenarios that we expect to broadly cover the phenomenology of the most well-motivated BSM models featuring new neutrino interactions below the MeV scale (these scenarios are summarized in Table I).
All scenarios correspond to different combinations of i) the number of interacting neutrino families, N int , ii) the internal degrees of freedom of the X particle, g X , and iii) if the X species has a non-zero primordial abundance or not, parametrized by ∆N BBN eff .To be specific, we consider the following: • Case (a), with N int = 3 and g X = 1, corresponds to the singlet majoron model in which neutrinos are pseudo-degenerate (note that pseudodegenerate neutrinos imply a universal coupling λ ν ).
• Case (b), with N int = 3 and g X = 3, corresponds to the commonly studied model of a light Z ′ boson coupled to a lepton number family symmetry.In this model it is once again a good approximation to consider a flavour universal coupling, since the PMNS matrix does not show a hierarchical structure.
Scenario Specification (a) • Case (d) corresponds to a case where a vector boson couples to a single neutrino mass eigenstate.As in scenario (c), this option is relevant in particular for 2m lightest ν < m X < 0.1 eV.However, a concrete model realization for m X > 0.1 eV in which a vector interacts only with one neutrino mass eigenstate is challenging, and generically involves cancellations of different couplings in flavour space.
• The cases (e) and (f) correspond to the cases (a) and (b), respectively, but allowing for a non-zero primordial abundance of the X particle parameterized by ∆N BBN eff .Such a primordial abundance of X particles can arise e.g.due to the decay of other, heavy particle species in the early Universe.For instance, majorons can be produced from the decays of GeV−scale sterile neutrinos [34], and the U (1) Lµ−Lτ gauge boson can be produced via muonantimuon annihilations in the early Universe [43].

III. COSMOLOGICAL IMPLICATIONS AND FORMALISM
Cosmological Implications: The cosmological implications of these light neutrino-philic bosons are governed by their decay rate into neutrinos.In particular, the ratio between the decay rate of X into neutrinos and the Hubble parameter at T ≃ m X /3 determines whether or not the X boson thermalizes in the early Universe.In a radiation dominated Universe, this ratio can be parametrized by: , where ⟨Γ(νν → X)⟩ is the thermally averaged inverse decay rate.For K eff ≳ 1 the X boson thermalizes with the neutrinos in the early Universe via decays and inverse decays out of neutrinos 2 .Thermalization has two important cosmological consequences: 1. Non-standard expansion at T ν ≲ m X -If the X boson thermalizes with neutrinos it will represent a non-negligible fraction of the energy density of the Universe.In particular, the X boson will behave as radiation until T ν ∼ m X but after it will start redshifting like matter and decay.This leads to a non-standard expansion history during this time, and to an enhanced value of N eff at the time of recombination (provided that m X has decays before recombination).
2. Suppression of neutrino free streaming -The new interactions between neutrinos and the X particle tend to homogenize the neutrino fluid, suppressing neutrino free streaming.This has important consequences for CMB observations as highlighted in the introduction.
Background Thermodynamics: The exact description of the thermodynamic evolution of the Universe in the presence of a light boson interacting with neutrinos can be found by solving the Liouville equation for the distribution function of neutrinos and the X boson.This is numerically very costly, but Ref. [49] explicitly demonstrated that for scenarios where the X boson interacts efficiently with neutrinos, namely for K eff ≳ 10 −3 , the thermodynamics can be accurately described by simple ordinary differential equations tracking the temperature and chemical potential of the neutrinos and the new light boson.These equations are explicitly outlined in Appendix I.In the left panel of Figure 2 we highlight the thermally averaged inverse decay rate (⟨Γ νν→X ⟩ = δρ X /δt| νν→X /ρ ν ) normalized to the Hubble parameter for a m X = 1 eV boson with g X = 1.We show the evolution for several values of K eff = 100 , 1 , 10 −2 representing cases where thermal equilibrium is well established, where thermal equilibrium is only slightly reached, and where the X boson does not thermalize, respectively.The energy density evolution for the X particle for each 2 Processes such as XX ↔ νν are only effective for λν ≳ 10 −7 and as can be seen from Eq. ( 5) we will be interested in much smaller couplings.of these cases is highlighted in the right panel of Figure 2.

Model
From this figure we can clearly see that for K eff ≳ 1 the X boson thermalizes with neutrinos and its thermodynamic evolution is dictated by thermal equilibrium.On the other hand, for K eff < 1 thermal equilibrium is not established which leads to out of equilibrium decays.The evolution at T ν ≲ m X /3 will lead in all cases to a nonstandard expansion history.
For K eff ≫ 1 and for m X ≳ 10 eV thermal equilibrium dictates what is the value of the neutrino energy density after the X particle has decayed away.By assuming thermal equilibrium and tracking the number and entropy densities of the neutrinos and X species (see [49]), we can calculate the minimum values of ∆N eff at recombination for the scenarios (a)-(d).These results are outlined in Table II.For X being a scalar mediator one expects ∆N CMB eff = 0.08 − 0.12 and for the vector mediator case ∆N CMB eff = 0.15 − 0.24.We note that these values are similar to Planck's 1σ sensitivity to N eff , and thus an accurate treatment of this modified expansion history is needed to analyze the latest data.
In the event that a primordial population of bosons already exists at the time of BBN, the process of thermalization at late times, i.e. near recombination, can significantly increase ∆N eff .For this reason, we differentiate the abundance of the new bosonic species at BBN and recombination using ∆N BBN eff and ∆N CMB eff .We illustrate the evolution of ∆N CMB eff assuming a primordial abundance of ∆N BBN eff = 0.4 in Figure 3. Two immediate conclusions can be drawn from this figure.Firstly, the shift in ∆N eff between BBN and recombination can greatly exceed the values outlined in Table II.Secondly, ∆N eff increases dramatically for K eff ≲ 1.This is because the X boson becomes non-relativistic and its delayed decay leads to a significant increase of the relative energy stored in this species.Consequently, scenarios with λ ν → 0 and ∆N BBN eff ̸ = 0 lead to a drastically distinct phenomenology compared to ΛCDM.Although, the effect of neutrinofree streaming suppression is negligible, these scenarios will be tightly constrained from the increase in ∆N eff .Cosmological Perturbations: In order to track the cosmological perturbations of the fluids describing neutrinos and the neutrino-philic boson X, we rely on several approximations.First, we treat the two interacting fluids as coupled, as done in past literature [30,31].This implies that we can evolve the perturbations jointly.In Left: Effective interaction rates at the background level (solid lines) as well as at the perturbation level (dashed lines) for different values of K eff .The scenario considered consists of all 3 neutrinos interacting with the scalar type boson.Right: Evolution of the normalized X boson energy density for the same scenarios as before.For reference, we highlight in dashed the photon energy density.
Evolution of N eff for the case of a scalar interacting with three neutrinos with a primordial contribution to ∆N BBN eff = 0.4.We notice that the value of N eff always increases and that for small K eff it increases significantly due to very out of equilibrium decays of the X particle.
the limit that the interactions are sufficiently strong this approximation is by definition valid.On the other hand, in the weak interaction limit, we also expect the approximation to be valid, because the perturbation equations in this case are equivalent to two decoupled fluids.
The second approximation adopted here enters in the collision term describing the 1 ↔ 2 interactions between the neutrinos and the X boson.Following Ref. [31] we assume: (1) Maxwell-Boltzmann statistics, (2) that the background momentum dependence of the neutrino distribution is not strongly time dependent, and (3) that the perturbation generated by gravity is universal to all the species involved.We expect all these approximations to hold in our scenario.
Finally, we treat neutrinos as being massless.This assumption significantly simplifies the evolution of the neutrino perturbations.Since current Planck data is con-sistent with massless neutrinos, setting an upper limit on the sum of neutrino masses at the level of m ν < 0.12 eV [4], we believe this approximation does not significantly alter our results.Nevertheless, a more thorough treatment including neutrino masses would be of interest, and thus we leave this for future work.
Under the approximations listed above, the equations describing the joint neutrino+boson system in synchronous gauge read [50]: Here, derivatives are taken with respect to conformal time, H is the conformal Hubble parameter, h and η represent the metric perturbations, a is the scale factor, ω = p/ρ is the equation of state of the system, c 2 s = dp/dρ is the sound speed squared, k defines the given Fourier mode, δ and θ are the energy and velocity perturbations respectively, F ℓ represents the ℓ moment of the perturbed distribution function, and the neutrino free streaming suppression rate is given by is [31]: In this expression we neglect the chemical potentials, which we explicitly checked to have negligible impact on observables.The coefficients are given by [31] where Γ(0, x) is the incomplete gamma function.At high temperatures Γ NF ∼ (m X /T ν ) 5 Γ(X → νν) and at very small temperatures Γ NF ∼ e −m X /Tν Γ(X → νν).This neutrino free streaming rate is shown as a function of temperature in dashed lines in the the left panel of Figure 2. We can clearly see that at high temperatures the scaling of Γ NF is different to the background evolution.Moreover at T ν ∼ m X /3, where the rate is maximal, it is a factor of ∼ 1/10 smaller than the background equivalent.It is actually easy to see that for Γ NF /H > 1, F ℓ → 0 exponentially fast, which strongly reduces neutrino free streaming.
Numerical Implementation in CLASS: We track the impact of the neutrino-X interactions on the CMB power spectrum by modifying the cosmological Boltzmann code CLASS [46,47].The code is available on github .It can also help to study the thermodynamic evolution of different BSM scenarios.
In the left panel of Figure 4 we show the evolution of the neutrino anisotropic stress associated with a mode of k = 0.1 Mpc −1 as a function of redshift.We choose k = 0.1 Mpc −1 because it is the largest wave number well probed by CMB observations.The evolution for different, smaller wave numbers are shown in Figure S8 of the appendix.From Figure 4 we can clearly see how the decays and inverse decays of X reduce the neutrino anisotropic stress.In the right panel of the same figure we also show the relative impact on the temperature power spectrum C T T ℓ compared to ΛCDM.The impact on the observable C T T ℓ spectrum can go well above the level of the 1σ relative error bars, as indicated by the grey band.
In Figure 5 we show the CMB temperature power spectrum for different values of m X , taking N int = 3, g X = 1, and fixing K eff = 10 4 .This corresponds to a scenario where the X particle interacts very efficiently with neutrinos, and thermal equilibrium is reached at T ∼ 30 × m X .From this plot we can appreciate a number of interesting features: firstly, we notice that for m X ≲ 0.1eV the impact on the CMB power spectrum is not significant.This is because the non-standard expansion history occurs after recombination, and owing to the high temperature suppression in the collision term, neutrino free streaming is not significantly altered before recombination.We notice that the most significant effect is for bosons with 1 eV ≲ m X ≲ 100 eV.This is because the interaction rate of these bosons is maximal during the window of redshift to which the CMB is sensitive, i.e. 2000 ≲ z ≲ 10 5 .Finally, for the case with heavy mediator, m X = 10 keV, the boson can not alter late-time free streaming, since it will have decayed already at higher redshift.This means that the observed effect purely corresponds to a shift in N eff of 0.12 (see Table II).

IV. CMB DATA ANALYSIS AND RESULTS
Cosmological Data and Analysis: We perform MCMC analyses with MontePython [51,52] on each of the models listed in Table I.For the likelihood we use data from Planck2018+BAO data [4,53].In particular, this includes the temperature and polarization power spectra, as well as the lensing likelihood, from Planck [53], and the 6DF galaxy survey [54], the MGS galaxy sample of SDSS [55], and the CMASS and LOWZ galaxy samples of BOSS DR12 [56][57][58][59].In order to investigate the extent to which these scenarios could explain or ameliorate the Hubble tension we perform additional MCMC analyses including a Gaussian likelihood on H 0 = 73.30± 1.04 km/s/Mpc [60].These results are used to replicate the three statistical criteria (described in detail below) introduced in the 'H 0 Olympics' [44].This comparison allows to establish the relative success and failure of the models of Table I in relation to other proposed solutions.
For the standard cosmological parameters and the nuisance parameters of the Planck likelihood we use the same priors as the Planck collaboration.For the mass and coupling of the neutrino-philic bosons we adopt log priors over the range: The lower bound on m X corresponds to twice the minimum mass of the heaviest neutrino, 2 |∆m 2 atm | ≃ 0.1 eV.For the case of the X boson interacting with N int < 3 neutrino families, the prior range is extended to log 10 (m X /eV) ∈ [−4, 3.5] as one of the neutrinos could be much lighter and thus open up parameter space for lighter X bosons.The lower limit in this case is chosen to be sufficiently small such that the interaction rate is never effective to thermalize the X boson.We also introduce a specific upper limit on λ ν = 10 −6 .This is because at larger couplings two-to-two processes (XX ↔ ν ν), which are not captured by our treatment, begin to become relevant.On the other hand, the lower limit in the coupling is chosen to be sufficiently small that the X boson is effectively fully decoupled from the neutrino sector.In this limit, λ ν → 0, ΛCDM is recovered.At sufficiently large masses, the X boson decays at high redshift, producing a shift in ∆N eff without altering neutrino free streamingour upper bound on the mass is set by the fact that this effect is the same for m X ≳ 1 keV (assuming a sufficiently large coupling such that the bosons thermalize).Finally, in some of the scenarios we also allow for a non-zero initial abundance of the X particle.We parameterize it by ∆N BBN eff and adopt a flat, linear prior over the range Performing the MCMC analysis with the likelihoods and priors as described above leads to the result of   5. Fractional difference on the TT power spectrum with respect to ΛCDM for the case of a scalar particle interacting efficiently with neutrinos, K eff = 10 4 , see Eq. ( 5).We show the results for different values of mX .
runs contain a total of N ∼ 2 × 10 6 samples.The 3σ exclusion region is obtained by binning the points in log 10 (m X /eV), and in each bin determining the coupling λ ν for which 99.7% of the samples have λ ν ≤ λ limit .A particularly interesting result is obtained for the scenario (c), i.e. the scalar boson X which interacts with N int = 1 neutrino family.In this scenario, we find a slight statistical preference for non-zero neutrino interactions; we note, however, that the ΛCDM limit is also favored at the 1σ level, implying the statistical preference for this best-fit region is not remarkably significant.This region can be seen more clearly in Figure S10, where the MonteCarlo samples are explicitly shown.This best fit region of parameter space roughly corresponds to: namely, this preferred region of parameter space cor-responds to scenarios where the neutrino anisotropic stress starts to be damped right before recombination, 1100 ≲ z ≲ 3500.This is highlighted by the red region labelled 'best fit region' in Figure 1.
We note that we do not find such a preferred region of parameter space for scenario (d) with a gauge boson interacting with a single neutrino species.The suppression of neutrino free-streaming is very similar to the case of a scalar, and thus we attribute the lack of preference for parameter space to the fact that the vector boson leads to a substantially enhanced expansion history for which Planck is sensitive to, see the lower row of Table II.

Implications for the Hubble Tension
It has been shown in [33][34][35] that models with neutrino X-boson interactions can have the potential to significantly ameliorate the Hubble tension for two main reasons: 1) the X-neutrino interactions can lead to a nontrivial enhancement of the expansion history near recombination, 2) there exists a level of degeneracy between the impact of the damping of neutrino free streaming and an enhanced value of N eff which allows for additional radiation without spoiling the fit to the data from Planck.In particular, the detailed statistical analysis of the 'H 0 Olympics' [44] awarded the model with a silver medal.However, as mentioned above, the original implementation of this model relied on numerous approximations.For this reason, we revisit the three 'H 0 Olympics' criteria using the improved analysis developed here.These criteria include: 1.The Gaussian Tension, given by where the minimum χ 2 is evaluated using a likelihood that does (C + SH 0 ES) and does not contain (C) the SH 0 ES likelihood.

Akaike Information Criterium (AIC), given by
where M refers to the model under consideration and N corresponds to the number of free parameters of that model.Here, the χ 2 min values are obtained using a likelihood that includes the Gaussian contribution from SH 0 ES.
Each criteria is intended to address a slightly different question -we refer the interested reader to [44] for a broader overview of the benefits and drawbacks of each.The results of each model are summarized in Table III.There we also show for comparison the ΛCDM result and the simple scenario containing free streaming dark radiation as parameterized by ∆N eff .Interestingly, none of the models investigated show a significant reduction in the cosmological tension, with the most successful of them only reducing it to the 3.2σ level (in comparison with 4.5σ for ΛCDM).This result obtained here represents a degradation compared to what was found in previous works [33][34][35].The main reason for this deviation is due to the refined collision term included here, see Eq. (7), which reduces the damping of neutrino free streaming with respect to the approximation of [33][34][35] at T ≫ m X .In particular, the full collision term helps to break the partial degeneracy between the damping of the neutrino free streaming at high redshift and the enhancement of ∆N eff .

V. ADDITIONAL CONSTRAINTS
The models we have discussed in the main text are subject to additional constraints coming from other cosmological probes, emission from astrophysical objects, and laboratory searches.In this section we briefly highlight the origin of each constraint shown in Figure 1.
Laboratory Constraints: In the two benchmark particle physics models we consider, see Eqns. ( 1)-( 2), the coupling of the new boson to neutrinos is constrained by a different set of laboratory constraints.In the case of X being identified as a light scalar, its coupling to neutrinos can give rise to double beta decay along the emission of a scalar.The latest constraints on λ ν from the non-observation of such a process from the EXO-200 experiment reads: λ ν < 0.9 × 10 −5 [61].In the case of X being a light U (1) Lµ−Lτ gauge boson, we adopt a nominal value of kinetic mixing induced at 1-loop by muons and taus, ϵ ≃ −g µ−τ /70 [62].The presence of this mixing can in turn change the scattering rate of neutrinos and electrons, which has been precisely measured by Borexino [63].For m X ≲ MeV, the coupling is constrained to be g µ−τ < 4 × 10 −5 [64][65][66].Both the EXO-200 and Borexino bounds are shown in Figure 1.
Supernova Bounds: Despite being very weakly coupled, the neutrino-philic bosons considered in this work can be copiously produced in extreme astrophysical environments such as supernovae.If so, these particles can modify the energy and temporal distributions of the neutrino flux arriving on Earth.In particular, in the majoron model the neutrino coalescence νν → ϕ can produce a delayed high-energy neutrino signal [67][68][69][70].The nonobservation of such a signature in the measured neutrino flux from SN1987A [71][72][73] leads to the following constraint [67]: for 10 keV ≲ m X ≲ 1 MeV.
On the other hand, the high densities present at supernovae induce flavour and helicity dependent effective neutrino masses.Therefore, for masses m X ≲ 10 keV, the process ν → νX in kinematically allowed [74,75].Including these processes one finds constraints at the level of The SN1987A bound for a U (1) Lµ−Lτ gauge boson were derived in [43,76].The emission of gauge bosons of m Z ′ < MeV is dominated by semi-Compton processes µγ → µZ ′ and the constraint imposed by the observation of the SN1987A signal is at the level of g µ−τ ≲ 10 −9 [76].
Star Cooling: A light U (1) Lµ−Lτ gauge boson with the canonical kinetic mixing interacts with charged matter, and thus can be produced in stars.Should these particles be produced, they can free stream out of the star, carrying away a sizeable amount of energy.Consequently, strong constraints can be derived by requiring that the stellar cooling rate is not significantly altered.Recasting the limits derived in [77] (see also [78] and [79]) using the nominal kinetic mixing ϵ = −g µ−τ /70 yields the bound in Figure 1, labelled 'Stars'.
BBN Bounds: The production of new relativistic particles prior to BBN will enhance the value of ∆N eff .This modifies the expansion rate and in turn the prediction of the primordial element abundances.Current observations of the primordial abundances are consistent with ∆N eff ∼ 0. In particular, ∆N BBN eff ≤ 0.41 at 2σ [80,81], and thus large deviations from this can yield strong constraints on the interactions with new particles.
Limits were recently derived on the majoron by identifying the couplings for which νν → ϕ lead to a shift in ∆N eff at the level of 0.5 [33].Comparable constraints were derived on the µ − τ gauge boson from the production of a primordial population via µ + µ − → Z ′ γ processes [43].These constraints are shown in Figure 1 with the label 'BBN'.

CMB bounds on out of equilibrium decays:
The thermodynamic treatment of the neutrino-philic bosons used in this study is only capable of accounting for moderate departures of thermal equilibrium, namely for K eff ≳ 10 −3 [49].In the absence of a primordial abundance, the region of parameter space with K eff ≲ 10 −3 is irrelevant as K eff controls the production of X particles and for such small K eff the energy density of X particles is negligible.However, even a small primordial abundance in the weakly coupled limit can yield strong observable consequences.The reason is that the primordial species can become non-relativistic prior to matterradiation equality, dramatically increasing the relative energy density stored in this species before it undergoes an out-of-equilibrium decay into neutrinos.The detailed treatment of this scenario is rather intricate (see e.g.[82,83]), and a full parameter space exploration is still lacking.In order to illustrate where these constraints would lie, we assume a primordial abundance at BBN of ∆N eff | BBN = g X × 0.027 (corresponding to the minimal value predicted for a boson that was in thermal equilibrium at temperatures above the electroweak phase transition) and derive an approximate constraint by requiring that N eff < 4 at recombination.We did this by tracking the evolution of the X boson energy density allowing for out of equilibrium decays and neglecting inverse decays (which are highly inefficient in this region of parameter space).In Figure 1 this constraint is indicated by the pink region labelled 'out of equilibrium decay' (and would exclude couplings below this line).

VI. SUMMARY, CONCLUSIONS AND OUTLOOK
In this work, we have presented an improved treatment of the cosmological evolution of weakly coupled neutrinophilic bosons with masses in the O(eV) range.This work represents a significant improvement upon previously analyses [33,34], which focused exclusively on the singlet majoron model and relied on a number of simplified approximations.Specifically, in this manuscript we present three updates: 1. We have incorporated the thermodynamic evolution tracing the out-of-equilibrium thermalization of the neutrino-philic bosons directly in the Boltzmann solver CLASS.This allows for a more accurate and careful treatment of the neutrino-boson interactions across a wide array of parameter space.The developed code is made public on github .
2. We have incorporated a recently derived collision term [30,31], which captures the impact of these interactions on the damping of the neutrino anisotropic stress.
3. We generalize this analysis to include: interactions with one, two, or three neutrino species, and both vector and scalar bosons.Our fiducial limits are recasted in the terms of the singlet majoron model and the U (1) Lµ−Lτ gauge boson, but these limits can be easily interpreted in the context of many other neutrino-philic boson models.
As shown in Figure 1, the limits derived using a combination of CMB and BAO data provide the strongest constraints to date across a range of masses near the O(eV) scale.We have also revisited the extent to which neutrino-philic bosons can resolve the Hubble tension.We show that the improved collision term, which is strongly suppressed in comparison to the previous approximations at T ≫ m X , significantly degrades the extent to which neutrino-philic bosons can ameliorate the tension.
In the case of the majoron singlet model, there exists a slight preference in the data for non-zero majoronneutrino interactions (at the ∼ 1σ level).This region of parameter space is expected to be fully probed in the near future by LiteBIRD [84] thanks to a cosmic variance limited measurement of the large scale EE polarization power spectrum.Upcoming observations from the Simons Observatory [85] are expected to measure N eff with a 1σ precision of 0.05.This will be an improvement by a factor of 4 as compared with Planck and will significantly improve sensitivity for bosons with masses 1 eV ≲ m X ≲ 1 MeV that thermalize in the early Universe with neutrinos.Both of these experiments are fully funded and expected to probe these regions of parameter space within a decade.
We use the formalism developed in [49,86] to trace the evolution of the background, which assumes that the distribution functions for all relevant species can be characterized by their temperature T i and chemical potential µ i .The time evolution equations for these quantities reads where in these expressions ρ i , n i , and p i are the energy density, number density, and pressure of the given species i.
The change in number and energy density of X bosons per unit time, δn X /δt and δρ X /δt, are given by Here we have used the Maxwell-Boltzmann approximation, and introduced the Bessel functions K i .Since the process we are considering is 1 → 2, the transfer rate for neutrinos are related to these shown here via: δρ ν /δt = −δρ X /δt and δn ν /δt = −2δn X /δt.For photons and neutrinos which are not coupled to the light neutrino-philic boson, one simply has dT dt = −H T .The evolution equations shown above, along with the modified perturbation equations outlined in Eqns.(6a)-(6d), are implemented into the publicly available cosmological Boltzmann solver CLASS [47].The incorporation of the new interactions into CLASS follows the standard methodology and is done in three steps.First, the input module is modified to read the new model parameters, including the coupling constant λ ν , the mass of the new boson m X , and its primordial abundance parametrized via ∆N BBN eff .A list with a detailed description of all readable parameters can be found in the file majoron.ini in the github repository of our code.The background evolution of the system as governed by Eq. (S1) is implemented into the function int background derivs within the background module.Its initial conditions are defined in the function int background initial conditions.Because the differential equations governing the background evolution are generally stiff, we explicitly implement a reduction of the step size when the energy density in the new species is non-negligble -specifically, we take this range to be 0.65 < T ν /m X < 10.The step size can be controlled in the .inifile via the parameter fine steps maj and is typically O(10 −5 ).The evolution of all relevant background quantities can be accessed by different modules through the pointer pba→.This is of particular importance to calculate the perturbations in the perturbations module.By accessing the background evolution, it is then straightforward to implement Eqns.(6a) -(6d) into the function int perturb derivs.

II. IMPACT ON COSMOLOGY
Here, we take the opportunity to provide a more detailed picture of how neutrino-philic bosons alter the temperature and polarization power spectrum.We begin by showing the evolution of the energy density stored in both neutrinos and the X boson for scenarios with different number of interacting neutrino species, N int , scalar or vector X-bosons, and with and without primordial X abundance.We then discuss the impact on the cosmological perturbations.In particular, we show the evolution of the speed of sound c 2 s , equation of state ω, the neutrino density contrast in the synchronous gauge δ ν , and the neutrino anisotropic stress σ ν , and finally the impact on the temperature and polarization power spectra C T T ℓ and C EE ℓ .In figure S6 we show the evolution of the energy density of the neutrino and X boson for various boson masses and interaction strengths, varying the primordial abundances, the number of interacting neutrinos, and the spin of the neutrino-philic boson.The general shape of the energy density evolution of the X boson is the same as long as Evolution of the energy density of the neutrino-philic boson X as well as the neutrinos for a variation of different cases as specified in the plot labels.
∆N BBN eff = 0.This can be seen from the left column of Figure S6.The only difference between the cases with different number of interacting neutrino families and X being a scalar or vector boson is the peak density the X boson reaches (while keeping the other model parameters fixed).In general, the more neutrino families interact with the X particles, the more the X particle will be populated in the thermal plasma.The same holds for the vector versus scalar case, i.e. ρ vector X > ρ scalar X at the peak of its thermalization history.Additionally, in the bottom right panel of Figure S6 we show the evolution of the X particle density for fixed K eff = 1 with varying mass m X .We see that max(ρ X ) is approximately the same for all depicted masses.On the other hand, if e.g.∆N BBN eff = 0.1, the relative change of the energy density between its initial value and its maximal value is relatively small.However, the striking feature is that for K eff ≪ 1 the boson energy density, ρ X , reaches larger values at lower temperature compared to the K eff ≫ 1 scenario.This is contrary to the case of vanishing primordial X particle abundance.The reason is, as outlined in Section III, that the X boson becomes non-relativistic and its small coupling leads to a delayed out-of-equilibrium decay.Lastly, in the top right panel of Figure S6 we also show a exemplary evolution of the neutrino energy density ρ ν for the case of X being of scalar type with vanishing primordial abundance and interacting with all three neutrino families.We choose to fix m X = 1 eV and vary K eff .This makes it evident that the energy density stored in the neutrinos is enhanced for K eff ≳ 1 compared to weakly and non-interacting neutrino scenarios.On the other hand, ρ ν (T ν → 0) also saturates to a maximal value and becomes independent of K eff as long as K eff ≳ 1.We can also translate the evolution of these energy densities into N eff , see Figure 3, via where ρ rad is the total energy density stored in radiation.This directly relates the late time enhancement of the energy density in the neutrinos and X boson to the observable measured by Planck.
The solutions derived from solving Eq.S1 also allow to compute the sound speed, c 2 s , and the equation of state, ω, of the joint fluid.The solution for two different masses of the neutrino-philic boson and different interactions strengths are shown in figure S7.We see that the maximal deviation from the relativistic approximation c 2 s = ω = 1/3  is reached for K eff = 1, independent of the mass of the X particle m X .In particular, we can appreciate that also for m X ≪ 1 eV the evolution significantly deviates from the relativistic approximation.Although differences of O(4%) can arise with respect to the relativistic approximation, we explicitly checked that this leads to a negligible effect on all observables.This approximation typically induces an error in the TT power spectrum at the level of O(0.01%) in all relevant regions of parameter space and is therefore well below Planck sensitivity.
We also show the evolution of the density contrast δ ν and the neutrino anisotropic stress (shear) σ ν for two different wavelengths in Figure S8.For reference, the black dotted lines indicate the ΛCDM expectation.The neutrino interactions lead to a significant reduction of the neutrino anisotropic stress at the time of recombination z ∼ 10 3 , see also Figure 4. On the other hand, the very same interaction leads to a significant enhancement of the density contrast at the same time.
Finally, we show the variation of the temperature fluctuation C T T ℓ and polarization C EE ℓ spectrum in the left and right panel of Figure S9.The temperature polarization spectrum shown here is complementary to the one in Figure 5 of the main text, in which we show the variation for fixed interactions strength and different masses m X .Here, we fix the mass to be m X = 1 eV and vary the interactions strength K eff .As expected, interactions with strength K eff ≫ 1 Fractional difference on the TT (EE) power spectrum with respect to ΛCDM for the case of a scalar particle interacting with neutrinos in the left (right) panel.We show the results for fixed mass mX = 1 eV and different values of K eff , see Eq. ( 5).lead to a significant perturbation of both spectra, C T T ℓ and C EE ℓ , which can exceed the 1σ error bars of the Planck mission shown in grey.In particular, the interactions induce a periodic perturbation spectrum with strong damping for the high-ℓ multipole moments as dictated by Eq. (7).
Having implemented the background and perturbation differential equations for the joined neutrino-X fluid, see Section I for details, an accurate parameter space scan can be done via a MCMC analysis.We chose to implement our code into the publicly available software MontePython [51,52].The full analysis then leads to the main result as shown in Figure 1.Here, we would like to give more specific details on two important aspects -i) how we have identified the exclusion region and the best-fit region, and ii) the correlation of neutrino interactions with different FIG.S10.Results of a MCMC analysis against the Planck legacy data.Black scattered points correspond N ∼ 2 × 10 6 Monte Carlo samples in the analysis.We find that 99.7% of all points are below the blue line, labeled as 3σ exclusion.In the left (right) panel we show the results for a scalar particle interacting with Nint = 1 (Nint = 3) neutrino families.Interestingly, for the Nint = 1 case the MCMC analysis identifies a 1σ preferred region, as can be seen by the red region.This corresponds to a non-trivial clustering of points which furthermore exhibit significant neutrino-philic interactions.In addition, for reference, we show in purple dashed the isocontours of fixed K eff .cosmological parameters.Let us start with point i).After removing the non-markovian points as well as the burn-in points of each chain, the raw points projected onto the parameter space of (m X , λ ν ) can be visualized as in Figure S10.For clarity we choose to only depict the case of X being a scalar mediator, but the same arguments also hold for the vector scenario.The exclusion region is simply found by demanding that within a given mass bin 99.7% of all points are below a given value of λ ν,i .We chose a bin size of ∆ log 10 (m X /eV) = 0.2.On the other hand, the best fit region can be obtained by evaluating the cluster density of sampled points.The same preferred region can also be found by running MCMC analysis softwares as e.g.GetDist [87].This brings us to point ii) -the triangle plots of different MCMC runs.We analyze the chains with the publicly available GetDist software.Of all the MCMC analysis done, we select a representative set of results.
• Comparison between the case of X being of scalar or vector type with fixed ∆N BBN eff = 0 and N int = 3.The likelihood to be tested against is the full Planck + BAO set.The result is shown in Figure S11.The lower bound on the mass of m X is set by the requirement of m X > 2 × m ν,i .We can see that for the X vector case globally slightly lower couplings λ ν are allowed compared to the X scalar case.Interestingly, the analysis indicates that the X vector case is compatible with slightly larger values of H 0 and is compatible with H 0 > 70 at the 2σ level for the region of high masses and low coupling.This is due to the contribution of the X particle to N eff as explained above and in the main text.
• The scenario of X being of scalar boson with ∆N BBN eff = 0 and N int = 1.The likelihoods are the same as before.The result of the analysis is shown in Figure S12.It clearly highlights the non-trivial 1σ preferred region in the plane (m X , λ ν ).In particular, this region indicates a slight preference for neutrino-X interactions such that the X boson starts reducing neutrino free-streaming by redshift z ∼ 1000 − 3500, see Eq. ( 13).
• The result for the scenario which allows a primordial abundance of the X scalar boson and N int = 1 is shown in Figure S13.Different combinations of ∆N BBN eff ≥ 0, g X = (1, 3) and N int = (1, 2, 3) lead to qualitative same results with only slight quantitative differences, as can be seen from Figure S6.The likelihoods are the same as before but now include also a) the Pantheon data set and b) the Pantheon data set together with the SH 0 ES prior.Contrary to what was found in Ref. [33], our refined analysis shows that even in the case b), see Table I, the H 0 value predicted by the model can not be increase to the 1σ SH 0 ES measured value.We find no significant increment in the prediction of the H 0 parameter.More details on the quantification of the H 0 -tension can be found in Section IV of the main text.
• For the case of a scalar interacting with N int = 3 neutrinos and with ∆N BBN eff = 0, we show in Figure S14 the full correlation of the standard cosmological parameters.The result is compared to ΛCDM in the same figure.We see that both cosmologies lead to similar correlations in these parameters modulo a small shift on H 0 and a multimodal posterior in the H 0 -ω CDM plane.

FIG. 4 .
FIG.4.Left panel: Evolution of the neutrino anisotropic stress for a mode of k = 0.1 Mpc −1 for ΛCDM and an scenario with Nint = 3 neutrinos interacting with a scalar with different coupling strengths.Right panel: Relative difference of the TT power spectrum in a majoron cosmology with respect to ΛCDM as a function of multipole ℓ.We show for reference the size of the Planck error bars.The comparison has been made with fixed standard cosmological parameters.We can clearly appreciate how the strong damping of the neutrino anisitropic stress on the left hand side is strongly related with a strong change on the power spectra.

2 eV c 2 s ω 2 FIG
FIG. S7.Speed of sound and equation of state of the joint neutrino+X boson system.A fully relativistic fluid will have c 2 s = ω = 1/3.

5 ×
FIG.S8.Density contrast δν and the neutrino anisotropic stress (shear) σν in synchronous gauge for fixed mass mX = 1 eV, but different wavelengths k and interaction strengths λν .

X 0 X 0 X
FIG. S11.1σ and 2σ posterior probabilities for the scenario of Nint = 3, ∆N BBN eff = 0.The case of X being a scalar boson is shown in blue and the vector case is shown in red.The model is tested against the likelihood of the full Planck18+BAO data set.
FIG. S14.1σ and 2σ posterior probabilities for the scenario of Nint = 3, ∆N BBN eff = 0 and X being a scalar boson.We show the correlation of the standard cosmological parameters and compare to ΛCDM.

TABLE II .
Minimum contributions to ∆N eff at the time of recombination resulting from the thermalization and subsequent decay of the X neutrino-philic boson.This corresponds to K eff ≫ 1 and mX ≳ 10 eV.
H 0i and σ i are the central value and the uncertainty on the inferred value H 0 .The index i = {C, SH 0 ES} refers to the cosmologically inferred value (using Planck and BAO) or the value measured by SH 0 ES, H 0 = 73.3± 1.04 km/s/Mpc. where