Probing the neutrino chemical potential with cosmological observations

The neutrino chemical potential has two main effects on cosmology, specifically on BBN and on the neutrino energy density, ofa different strength depending on the neutrino flavour. The BBN effect is mostly due to the fact that the electron neutrino chemical potential, $\xi_{\nu_{\mathrm{e}}}$, modifies the neutron-to-proton conversion rates during BBN [19], shifting the neutron–proton equilibrium ratio at freeze-out and thereby increasing (decreasing), for negative (positive) $\xi_{\nu_{\mathrm{e}}}$, the primordial Helium-4 abundance, expressed by the fraction $Y_{\mathrm{He}}$ [20]

Moreover, all chemical potentials $\xi_{\nu}$ enter the distribution function of cosmological neutrinos, expressed (for neutrinos and antineutrinos) as

where $y=q/T_{\nu}$ is the dimensionless comoving momentum. This effect changes the neutrino energy density, modifying the background expansion history at all times. In the early Universe this affects the effective number of neutrinos $N_{\mathrm{eff}}$,

At late times, the change in $\Omega_{\nu}$, determined by the change in the neutrino number density $n_{\nu}$, reads as

To illustrate the effects of a time varying chemical potential on the CMB and on the matter power spectrum, we refer to Figure˜2. We focus on three degenerate massive neutrinos and explore four models with time-varying chemical potential. Specifically, at four nodes in redshift (named $\xi_{\nu,\,\mathrm{late}},\,\xi_{\nu,\,115},\,\xi_{\nu,\,\mathrm{rec}},\,\xi_{\nu,\,\mathrm{BBN}}$, located respectively at $z=10,\,115,\,1100$ and $3\times 10^{8}$; see Table˜1 in Section˜III.1), $\xi_{\nu}(z)$ transits from 0 to 1, and then to 0 again. We compare these cases to those in which the chemical potential is constant, either 0 or 1, representing the values $\xi_{\nu}$ can assume in our (physical) limiting cases.

As we are not compensating any of the effects that an increase of $\xi_{\nu}$, translated in an increase of $N_{\mathrm{eff}}$ or $\Omega_{\nu}$, have on the expansion history, the models $\xi_{\nu,\,\mathrm{rec}}$ and $\xi_{\nu}=1$ (purple and yellow curves, respectively) undergo matter-radiation equality at redshift $z_{\mathrm{eq}}$ around $20\%$ smaller with respect to the rest of the models, due to the change to the energy density in relativistic matter around recombination. The change in the evolution of the neutrino energy density due to a time-depending chemical potential can be seen in the top-right panel of Figure˜2. Its effects on CMB $C_{\ell}^{TT}$ linear temperature power spectrum and the linear matter power spectrum $P(k)$ are, instead, displayed in the lower panels of the same Figure. Therein, we can appreciate how changes at small scales (large $\ell$ or $k$) are present if $\xi_{\nu}(z)$ deviates from 0 at early times, while if $\xi_{\nu}$ differs from 0 at late times, effects on large scales are present.

Matter power spectrum $\,-\,$ In the bottom-right panel of Figure˜2 we show the impact of a time-varying neutrino chemical potential on the matter power spectrum. In the $\xi_{\nu,\,\mathrm{late}}$ model (blue curve) we observe a suppression at large scales due to an increased Hubble parameter at late times, while on small scales, where free-streaming is present, it displays an increased step-like suppression due to an increased $\Omega_{\nu}$, an effect similar to the one in the $\xi_{\nu,\,115}$ case (green curve). As anticipated, the models $\xi_{\nu,\,\mathrm{rec}}$ and $\xi_{\nu}=1$ undergo matter-radiation equality at a later time, and this pushes the peak of the matter power spectrum to larger scales, explaining the suppression at small scales. Finally, the effect of $\xi_{\nu,\,\mathrm{BBN}}$ (light-blue curve) on $P(k)$ amounts to an increase of $N_{\mathrm{eff}}$ in the early times while keeping other quantities (namely the matter radiation equality and the baryon and cold dark matter abundances) fixed. This translates into an increase of $P(k)$ for scales smaller than the matter-radiation equality one.

CMB temperature power spectrum $\,-\,$ Regarding the $C_{\ell}^{TT}$ (bottom-left panel of Figure˜2), the slight deviation of $\xi_{\nu,\,\mathrm{late}}$ is too small to be visible at this scale. The deviation of $\xi_{\nu,\,115}$, instead, is due to a larger early ISW with respect to $\xi_{\nu}=0$ (black curve); contrary to the standard late ISW, the effect is not a tangible enhancement of the first peak, but a slight enhancement of power on larger scales. The effect of $\xi_{\nu,\,\mathrm{rec}}$ is due to the different $z_{\mathrm{eq}}$ in the model, shifting the angular diameter distance at decoupling and thus changing the position of the first peak and the envelope of all the secondary peaks, analogously to what an increased $N_{\mathrm{eff}}$, without compensating anything else in the expansion history, would generate. The modifications generated by the transition of $\xi_{\nu,\,\mathrm{BBN}}$ are purely induced by the (anti-)correlation between the primordial Helium fraction $Y_{\mathrm{He}}$ and $\xi_{\nu}$. This is the only effect which is sensitive to the sign of the chemical potential. As a matter of fact, increasing $\xi_{\nu}$ at BBN implies reducing $Y_{\mathrm{He}}$, which increases the CMB temperature power in the damping tail as fewer baryons are bound in Helium, leaving more free electrons before recombination. This shortens the photon mean free path, leading to a reduction of the diffusion damping.

In both spectra, we can see a combination of all the effects described above for the intermediate cases for the $\xi_{\nu}=1$ model (yellow curve).

One of the best methods to perform an agnostic data-driven determination of a particular function is via the Piecewise Cubic Hermite Interpolating Polynomial (PCHIP) [24, 25], which can preserve the shape of data distribution due to its special mathematical properties and avoid many of the oscillations and overshoots that ordinary cubic splines can produce. Overall, this method consists of a combination of third order polynomials that ensures smoothness between a set of input nodes. For more details on the PCHIP interpolation method, see for example the recent works of Refs. [26, 27].

In this work, we apply this method to reconstruct the neutrino chemical potential, $\xi_{\nu}$, in a model-independent way, using four nodes at fixed redshifts. In this way, the functional form is completely determined once the values of the chemical potential at the four redshift values of our choice are specified.

The choice to work with only four nodes is dictated by the result of a previous simulation done with a larger (seven) set of nodes and by computational time: a higher number of nodes would translate into a higher number of parameters to be sampled in the MCMC. For the test simulation, we considered seven nodes at redshifts $z=0,\ 1,\ 10,\ 115,\ 1100,\ 3\times 10^{8},\ 10^{14}$. Given that data are not constraining the late time nodes (i.e. $z=0,\ 1,\ 10$), as the effects of the neutrino distribution function are weaker at such epochs, we decided to collapse ³³3A clarification is in order: with “collapse” we mean that we still have an interpolation on seven nodes but we impose that the variables of interest in the ”collapsed” nodes assume the same value. them into a single node, specifically considering the one at $z=10$. We also collapse the last two nodes, $z=3\times 10^{8},\ 10^{14}$ into a single one to avoid time-varying $\xi_{\nu}$ during BBN. In other words, in our PCHIP runs $\xi_{\nu}$ is constant below $z=10$ and above $z=3\times 10^{8}$, respectively with values $\xi_{\nu,\,\mathrm{late}}$ and $\xi_{\nu,\,\mathrm{BBN}}$.

Therefore, as previously anticipated, we keep the nodes at redshift $z=10,115,\,1100,\,3\times 10^{8}$ to take into account the epochs at which a non-zero chemical potential is expected to have a more significant impact. Respectively, these nodes represent the epoch in which structure formation becomes non-linear, the epoch in which a neutrino with mass of around $0.06$ eV transitions from being relativistic to non-relativistic, the epoch of recombination and the time of BBN. The full set of nodes considered is shown in Table˜1.

We shall always assume a spatially flat Universe with adiabatic initial conditions. We consider two different scenarios, regarding how we treat neutrinos:

Case $\bm{3\,\nu}$$\,-\,$ We consider three degenerate massive neutrinos, meaning that we are assuming $\xi_{\nu_{\mathrm{e}}}=\xi_{\nu_{\mu}}=\xi_{\nu_{\tau}}\equiv\xi_{\nu}$. We shall also fix the sum of the neutrino masses to be $\sum m_{\nu}=0.06$ eV, which is the lowest limit allowed by neutrino oscillation experiments [28, 29, 30].

Case $\bm{(1+2)\,\nu}$$\,-\,$ We consider a very light neutrino state (m${}_{\mathrm{lightest}}=10^{-5}$ eV), that can be connected to the electron neutrino, as it would be the case for the normal hierarchy, with chemical potential $\xi_{\nu_{\mathrm{e}}}$. We then consider the remaining two neutrinos to be degenerate, massive ($\sum m_{\nu,\,\mathrm{degenerate}}=0.06$ eV) and with the same chemical potential, $\xi_{\nu_{\mu}}=\xi_{\nu_{\tau}}\equiv\xi_{\nu_{\mathrm{x}}}$ ⁴⁴4The names of the two chemical potential are given consistently with the nomenclature used in PArthENoPE..

In principle, each neutrino family could have a different chemical potential. The effect of neutrino oscillations, however, tends to re-equilibrate the difference between the three neutrino families [31] even when the initial values are quite apart from one another. In this sense, it is perfectly reasonable to assume that at least two neutrino families share the same degeneracy parameter.

We introduce an additional classification, regarding how we treat the neutrino chemical potential $\xi_{\nu}$. Specifically, we consider the two following scenarios:

Case $\bm{\xi_{\nu,\,\mathrm{constant}}}$$\,-\,$ In this case, we simply consider the neutrino chemical potential to be constant over time, leaving it as a free parameter in the MCMC analyses.

Case $\bm{\xi_{\nu,\,\mathrm{PCHIP}}}$$\,-\,$ In this scenario, instead, we reconstruct the chemical potential using the PCHIP interpolation method with four nodes, as described in the previous section.⁵⁵5Note that in the $\bm{\xi_{(1+2)\nu}}$ case, we reconstruct both $\xi_{\nu_{\mathrm{e}}}$ and $\xi_{\nu_{\mathrm{x}}}$. The aim of this reconstruction is not to set bounds on $\xi_{\nu}$ at different epochs, but to prove that BBN is the most sensitive epoch to $\xi_{\nu}$, testing the origin of the bounds one usually derives on the chemical potential by performing a cosmological analysis.

In the following we shall consider all the possible combinations of the four scenarios described above, namely we consider the cases $\bm{\xi_{3\nu,\,\mathrm{constant}}}$, $\bm{\xi_{3\nu,\,\mathrm{PCHIP}}}$, $\bm{\xi_{(1+2)\nu,\,\mathrm{constant}}}$ and $\bm{\xi_{(1+2)\nu,\,\mathrm{PCHIP}}}$.

To constrain the parameters of our models, we use a modified version of the Cosmic Linear Anisotropy Solving System code (CLASS)⁶⁶6The modified version of CLASS that supports the findings of this article is not publicly available, but it can be obtained from the authors upon reasonable request. [32, 33]. The implemented modifications aim to compute the neutrino distribution function for the case of a non-zero chemical potential, to include the PCHIP reconstruction method and to modify the interpolation performed on the BBN tables. Specifically, the BBN tables currently in CLASS do not account for a non-zero neutrino chemical potential. To obtain the new table, we used the BBN code PArthENoPE [34, 35, 22]. We then performed the inference of the cosmological parameters using Cobaya [36, 37] with the convergence set as a Gelman-Rubin test [38] of $R-1.$ For the cases $\xi_{\nu,\mathrm{constant}}$ our aim is to place limits, and therefore we require a convergence of $R-1\leq 0.01$, while for the scenarios $\xi_{\nu,\mathrm{PCHIP}}$ where our aim is to understand the source of the constraining power, we relax the threshold at $R-1\leq 0.04$. The posterior distributions are obtained and analyzed using the GetDist package [39].

Finally, we perform a Bayesian model comparison. For each model $\mathcal{M}_{i}$, we compute the Bayesian evidence⁷⁷7We actually compute its average value, with the associated error, obtained considering different values for the parameters burnlen, controlling the burn-in region of the chains, and thinlen, controlling the thinning, thus the markovianity, of the chains. using the MCEvidence package [40, 41], defined as

Specifically, we make use of a Cobaya-friendly version of MCEvidence as presented in the wgcosmo Github repository⁸⁸8https://github.com/williamgiare/wgcosmo/tree/main/statistics/MCMC_Evidence. We compare the models considered in this work with a $\Lambda$CDM-like model, where we consider $\xi_{\nu}=0$ but assume massive neutrinos, distinguishing between the $3\nu$ and $(1+2)\nu$ cases.

Model comparison is performed computing the logarithm of the Bayes factor

A positive (negative) $\ln\mathcal{B}_{ij}$ means that our $\Lambda$CDM model $\mathcal{M}_{i}$ is favored (disfavored) over the specific model $\mathcal{M}_{j}$ considered in this work. Using the revised Jeffreys’ scale [42], the evidence is classified as inconclusive⁹⁹9This means that the two models perform comparably. ($0\leq|\ln\mathcal{B}_{ij}|<1$), weak ($1\leq|\ln\mathcal{B}_{ij}|<2.5$), moderate ($2.5\leq|\ln\mathcal{B}_{ij}|<5$), strong ($5\leq|\ln\mathcal{B}_{ij}|<10$) and very strong ($|\ln\mathcal{B}_{ij}|\geq 10$).

In all the analyses performed, we sample over the usual $\Lambda$CDM set of parameters, i.e. $\{\omega_{b},\,\omega_{c},\,100\,\theta_{s},\,\tau_{\rm reio},\,\ln(10^{10}A_{s}),\,n_{s}\}$, to which we add the additional parameters characteristic of each cosmological scenario considered, as explained in Section˜III.2 and shown in Table˜2, together with the priors imposed on all the parameters. Note that all analyses are conducted with a fixed $N_{\mathrm{eff}}^{\mathrm{std}}=3.044$, where we actually mean that we keep fixed some of the CLASS parameters governing the neutrino sector, specifically deg_ncdm, N_ncdm and N_ur. As can be seen from Equation˜5, the total $N_{\mathrm{eff}}$ will be given by $N^{\mathrm{tot}}_{\mathrm{eff}}=N_{\mathrm{eff}}^{\mathrm{std}}+\Delta N_{\mathrm{eff}}(\xi_{\nu})$ and, in the PCHIP case it will not be a constant during cosmic evolution since the contribution given by $\xi_{\nu}$ will vary over redshift.

We exploit the current state-of-the-art cosmological datasets to place constraints on the cosmological scenarios described above. In particular, we consider two different types of measurements as the baseline dataset:

• •

  CMB measurements: we use CMB temperature data from the Commander likelihood [1, 43] and $E$-mode polarization data from SRoll2 [44, 45] at large angular scales ($2\leq\ell\leq 29$), complemented at high multipoles by the Plik_lite likelihood from the Planck mission. These are combined with the ACT-lite¹⁰¹⁰10https://github.com/ACTCollaboration/DR6-ACT-lite likelihood for ACT DR6 [46] and the SPT-lite¹¹¹¹11https://github.com/SouthPoleTelescope/spt_candl_data likelihood for SPT-3G D1 [47, 48]. Following the prescriptions described in Refs. [46, 47], we restrict ourselves to Planck data up to multipoles $\ell_{\rm max}=1000$ in temperature and $\ell_{\rm max}=600$ in polarization, and assume no correlations between SPT and the other datasets. Furthermore, we include the joint Planck+ACT DR6 lensing likelihood¹²¹²12https://github.com/ACTCollaboration/act_dr6_lenslike [49, 50, 51, 52] and the MUSE lensing likelihood¹³¹³13https://github.com/qujia7/spt_act_likelihood for SPT-3G D1 [53]. In the following, we shall refer to this combined dataset combination as SPA.

• •

  BAO measurements: we include the three-year data collection of Baryon Acoustic Oscillations (BAO) measurements from the Dark Energy Spectroscopic Instrument (DESI) presented in the second data release (DR2) [54]. We refer to this dataset as DESI.

Using the combination SPA + DESI as a baseline allows us to get constraining power from the early to the late Universe.

Since a non-zero chemical potential in the early Universe will impact BBN, we include two additional likelihoods based on Deuterium and Helium-4 abundances. Specifically, we compare the predictions on the Deuterium abundance obtained from PArthENoPE to the PDG value [23], $D/H=(25.47\pm 0.29)\times 10^{-6}$. On the other hand, we separately consider two different measurements of the Helium fraction, $Y_{\mathrm{He}}$, as listed in the following:

• •

  EMPRESS: recent results by the EMPRESS survey [17] report a value of the Helium fraction, $Y_{\mathrm{He}}=0.2370^{+0.0034}_{-0.0033}$, which is $\sim 3\sigma$ lower than the standard BBN (sBBN) prediction, $Y_{\mathrm{He}}=0.24709\pm 0.00017$ [21]. As a consequence, a non-zero value for the electron neutrino chemical potential is derived by the collaboration, $\xi_{\nu_{\mathrm{e}}}=0.05^{+0.03}_{-0.02}$.

• •

  LBT: the LBT $Y_{p}$ project has recently reported a new measurement of the Helium fraction [55, 16, 56], which does not present a tension with the standard BBN predictions. The value measured by the LBT collaboration, $Y_{\mathrm{He}}=0.2458\pm 0.0013$, is the most precise determination of the primordial ⁴He mass fraction to date.

We also take into account the theoretical errors introduced by the computations of PArthENoPE, as reported in [22], meaning that as reference (theoretical) values in the likelihood we will have $Y_{\mathrm{He}}$ = value $\pm\,\sigma\,\pm\,0.00003\,\pm 0.00012$ and $D/H$ = (value $\pm\,\sigma\,\pm\,0.06\,\pm 0.03)\times 10^{-5}$.

Using the importance sampling method implemented in Cobaya, we re-evaluate our SPA + DESI chains to consider the constraining power of BBN on our models.

Figure˜3 shows the marginalized posterior distributions for a reduced set of parameters sampled in the MCMC, in addition to the neutrino chemical potential.

When BBN measurements are not included, the contours enlarge, giving more freedom for $\xi_{\nu}$ to deviate from zero. It is interesting to notice that the posterior of $\xi_{\nu}$ is asymmetric: data prefer positive values over negative ones. This is a direct consequence of the effect that $\xi_{\nu}$ induces on $Y_{\mathrm{He}}$ during BBN; CMB data (and BBN data from EMPRESS) have a preference for slightly lower $Y_{\mathrm{He}}$ with respect to the LBT data on BBN, pushing $\xi_{\nu}$ to explore a region away from zero, see Table˜3.

When BBN likelihoods are included in the analyses, whether based on EMPRESS or LBT measurements, the bounds tighten significantly and the limits on the Helium fraction converge toward the respective collaboration measurements. As expected, the tightest bounds are obtained when including the LBT measurement of the $Y_{\mathrm{He}}$, as it is currently the most precise determination of the primordial Helium fraction available.

It is interesting to note that, even though the EMPRESS measurement creates a mild tension with standard BBN predictions, the results obtained in that case are consistent with our baseline (SPA + DESI) scenario.

We have also computed the Bayesian evidence with respect to a $\Lambda$CDM model, see Table˜3. Our results show that there is only weak evidence for the former model with respect to models with non-vanishing lepton asymmetries. There is only a moderate evidence for a $\Lambda$CDM Universe once LBT data are included in the analyses, which makes sense, given the fact that this dataset does not present any tension with the standard BBN prediction, and we find that cosmology alone prefers lower values for the Helium fraction, $Y_{\mathrm{He}}$.

We finally highlight that, as expected, $N_{\mathrm{eff,\,BBN}}$ is a direct function of $\xi_{\nu}$, and follows the predicted parabola shape (see Equation˜5).

Note that we can only derive an upper bound on $N_{\mathrm{eff,\,BBN}}$ since we always work with a fixed value for $N_{\mathrm{eff}}^{\mathrm{std}}=3.044$, meaning that the impact of $\xi_{\nu}$ is only to increase its value¹⁵¹⁵15This comment remain true for all the considered analyses.. Moreover, if we focus on the constraints presented in Table˜3 for $\Delta N_{\mathrm{eff,\,BBN}}$, we find them to be tighter than the ones presented by, e.g., the Planck collaboration. This is due to the fact that in our analysis, both $N_{\mathrm{eff,\,BBN}}$ and $\Delta N_{\mathrm{eff,\,BBN}}$ are derived parameters from $\xi_{\nu}$ (see Equation˜5), characterized by a uniform prior. As a result, the induced prior on both parameters is non-uniform, assigning a larger prior volume to lower values of $\Delta N_{\mathrm{eff,\,BBN}}$.

In Figure˜4 we show the marginalized posteriors for the same set of parameters as in the previous subsection, but distinguishing now between the lightest neutrino contribution and that of the massive species, treated as degenerate. This leads to splitting the previous $\xi_{\nu}$ into two parameters, namely $\xi_{\nu_{\mathrm{e}}}$ and $\xi_{\nu_{\mathrm{x}}}$, following PArthENoPE nomenclature.

The most interesting result is the bimodal posterior for $\xi_{\nu_{\mathrm{x}}}$, leading to “banana-shaped” contours in all panels involving this parameter, which is a direct consequence of the behaviour shown in Figure˜1. Moreover, for the remaining of the effects on cosmology, the bimodal distribution is explained by the fact that in Equations˜5 and 6 the sign of each $\xi_{\nu_{\alpha}}$ is physically irrelevant. As a consequence, the data cannot distinguish between positive and negative values of $\xi_{\nu_{\mathrm{x}}}$, but it can for $\xi_{\nu_{\mathrm{e}}}$ due to its important role in the weak interaction rates important for BBN.

Comparing Figure˜4 to Figure˜3, it is clear that, in the degenerate scenario, the bounds on $\xi_{\nu}$ are actually driven by the bounds on $\xi_{\nu_{\mathrm{e}}}$, while almost the whole range is allowed for $\xi_{\nu_{\mathrm{x}}}$, meaning that BBN (and cosmology more broadly) is mostly sensitive to $\xi_{\nu_{\mathrm{e}}}$. This is expected, as the electron neutrino is directly involved in BBN processes, specifically shifting the neutron–proton equilibrium ratio at freeze-out and thereby modifying the value of the primordial Helium abundance, whereas the other flavours mostly affect cosmology through their relativistic and non-relativistic energy densities.

In comparison with the previous scenario (e.g., three degenerate massive neutrinos), there is no significant change regarding the impact of the different BBN likelihoods considered on the parameters, with a general tightening of the posteriors, more evident when we consider LBT. Moreover, the LBT likelihood has also the effect of suppressing the bimodal behavior of the $\xi_{\nu_{\mathrm{x}}}$ posterior, making it converge to a single peak centred in 0. This is due to the high precision of the LBT measurement and its agreement with the standard BBN value for the Helium fraction, effectively recovering the standard predictions and removing the degeneracies between the parameters. In contrast, the EMPRESS measurement prefers lower value of $Y_{\mathrm{He}}$, leaving more freedom for deviations of $\xi_{\nu_{\mathrm{x}}}$ from 0, thus not suppressing the bimodal distribution. It also allows more freedom to $\xi_{\nu_{\mathrm{e}}}$.

A notable feature of this scenario is the appearance of a tail in the posterior distribution of $H_{0}$, which relaxes its upper bound towards larger values, possibly alleviating the tension with Supernovae measurements [57]. We argue that this tail is due to the larger values of $N_{\mathrm{eff}}$, compared to standard predictions, driven by the regions away from zero that $\xi_{\nu_{\mathrm{e}}}$ and, most significantly, $\xi_{\nu_{\mathrm{x}}}$ explore. This reflects the positive correlation between $H_{0}$ and $N_{\mathrm{eff}}$: a larger expansion rate in the early Universe reduces the sound horizon and keeping the angular acoustic scale constant requires an increase in the value of $H_{0}$. Notice that, compared to the previous scenario, the correlation between $\xi_{\nu_{\mathrm{e}}}$ and $H_{0}$ is positive, implying therefore a larger value of $H_{0}$, as the value of $\xi_{\nu_{\mathrm{e}}}>0$ (see Table˜4).

We highlight that, in this scenario, $N_{\mathrm{eff,\,BBN}}$ is a function of both $\xi_{\nu_{\mathrm{x}}}$ and $\xi_{\nu_{\mathrm{e}}}$, and has a large tail towards higher values due to the wide region explored by $\xi_{\nu_{\mathrm{x}}}$. The correlation between the latter and $\xi_{\nu_{\mathrm{e}}}$, in which $\xi_{\nu_{\mathrm{e}}}$ is only varying within a small region is responsible for the narrowness of the parabola in the correlation with $N_{\mathrm{eff,\,BBN}}$.

Finally, considering the Bayesian evidence, we find a (non-signifcant) preference for our model over $\Lambda$CDM when we include the EMPRESS likelihood, otherwise $\Lambda$CDM is mildly preferred (see Table˜4).

In Figure˜5 we show the reconstruction of the neutrino chemical potential, obtained with the PCHIP method. We only show the evolution between the nodes listed in Table˜1, since earlier than $z=3\times 10^{-8}$ and later than $z=10$ the chemical potential has the same values as the respective limiting nodes (see Section˜III.1). We show both the 68% and 95% C.L. results, as dark-blue bands for our baseline data set and lines for the post-processing cases, respectively light-blue solid for including EMPRESS and cyan dashed for including LBT. In the top panel we show the reconstruction for the common $\xi_{\nu}$ of the $3\nu$ case, while in the middle and bottom panel we show the reconstructions for the two chemical potentials introduced in the $(1+2)\nu$ case, respectively $\xi_{\nu_{\mathrm{e}}}$ and $\xi_{\nu_{\mathrm{x}}}$.

Case $\bm{3\,\nu}$$\,-\,$As one would expect, we find that the most constraining epoch is around the BBN period (see Figure˜5). This is due to the correlation between the chemical potential and the Helium fraction, to which CMB data are sensitive through the diffusion damping of photons, mostly affecting BBN predictions (see Figure˜1). Adding BBN data makes the constraint significantly tighter at the BBN epoch, propagating the effect up to today, still leaving more freedom to the values that $\xi_{\nu}$ can take.

Interestingly, we find that the chemical potential is constrained, albeit not tightly, around recombination, when neutrinos hold about one-tenth of the total energy density of the Universe. We argue that this constraint derives from the massive nature of neutrinos, and the fact that variations in the chemical potential reflect onto $\Omega_{\nu}$, to which CMB spectra, in particular CMB lensing, are sensitive through e.g. the free-streaming. This propagates a bit at later times, as previously stated, albeit without adding a significant constraining power to set tight bounds on the parameter $\xi_{\nu}$.

In Figure˜6 we show the 68% and 95% C.L. triangle plot on the nodes used for the reconstruction and, also, on their correlations with $H_{0}$, $Y_{\mathrm{He}}$ and $N_{\mathrm{eff,\,BBN}}$.

Focusing on the nodes, we see the emergence of bimodal distributions for $\xi_{\nu,\,\mathrm{rec}}$ and $\xi_{\nu,\,115}$, in a similar fashion to the bimodal distribution $\xi_{\nu_{\mathrm{x}}}$ (see Figure˜3). The reason for the bimodal posterior is an effect due to a mild preference for a larger $N_{\mathrm{eff}}$ around recombination (see also Ref. [58]), and for a slightly larger $\Omega_{\nu}$: such preferences are translated into values of $\xi_{\nu,\mathrm{rec}}$ and $\xi_{\nu,\mathrm{115}}$ away from 0, see Table˜3. The inclusion of the BBN measurements has the strongest impact on the BBN node, propagating lightly to the other nodes by killing the bimodal posterior distributions (particularly if we consider the LBT data, which drive $\xi$ to 0), while having zero effect on the late time node, consistently with our expectations.

Since we have more degeneracies between parameters and a time varying chemical potential, and therefore a time-varying $N_{\mathrm{eff}}$, we find a longer tail in the posterior distribution of $H_{0}$, partly lost once we include the BBN data.

Finally, the results for the Bayesian evidence with respect to a $\Lambda$CDM model are shown in Table˜3. In this case, our results show that there is only weak evidence for the PCHIP lepton asymmetry model with respect to the $\Lambda$CDM scenario when considering the basic dataset and also when adding EMPRESS observations to the baseline data. As in the constant case, there is only a mild moderate evidence for a $\Lambda$CDM Universe once LBT measurements are included in the analyses.

Case $\bm{(1+2)\,\nu}$$\,-\,$

It is particularly interesting how the effects driving the bounds - BBN and $N_{\mathrm{eff}}$ in the early Universe, free-streaming through $\Omega_{\nu}$ in the late Universe - decouple, once we allow for a separation between $\xi_{\nu_{\mathrm{e}}}$ and $\xi_{\nu_{\mathrm{x}}}$. During BBN, $\xi_{\nu_{\mathrm{e}}}$ directly modifies the neutron-to-proton conversion rates, and is by far the most constrained quantity, with similar behaviour as the $\xi_{\nu}$ of the previous case, while $\xi_{\nu_{\mathrm{x}}}$, which only enters BBN indirectly through $N_{\mathrm{eff}}$, is much less constrained by BBN light element abundances. Both quantities reach similar bounds around recombination, when the mass is still negligible. At later times, the very light neutrino state and the corresponding $\xi_{\nu_{\mathrm{e}}}$ which in practice does not contribute to $\Omega_{\nu}$, is much more loosely constrained with respect to $\xi_{\nu_{\mathrm{x}}}$.

Adding BBN datasets, as one expects, modifies the profile of the $\xi_{\nu_{\mathrm{e}}}$ around BBN, pushing it towards slightly positive values, and the propagation of this modification to later epochs and to $\xi_{\nu_{\mathrm{x}}}$ is small. As a confirmation that the data is only sensitive to $\xi_{\nu_{\mathrm{x}}}$ through even powers, the contours of $\xi_{\nu_{\mathrm{x}},\,i}$ are symmetrical at all epochs.

In Figure˜7 we show the 68% and 95% C.L. triangle plot on the nodes used for the reconstruction and, also, on their correlations with $H_{0}$, $Y_{\mathrm{He}}$ and $N_{\mathrm{eff,\,BBN}}$.

Focusing on the nodes, we see the emergence of bimodal distributions for the early $\xi_{\nu_{\mathrm{x}}}$ nodes, given the preference for larger integrated values, which are sensitive to even powers of the chemical potentials. However, the bimodal distribution does not appear in the first nodes of $\xi_{\nu_{\mathrm{e}}}$, due to its asymmetric effect in BBN physics. Adding BBN measurements has the same effect as in the previous case: it affects mostly the BBN epoch nodes, pushing $\xi_{\nu_{e},\mathrm{BBN}}$ closer to the values measured by each collaboration. Curiously, before adding these measurements, the posterior for $\xi_{\nu_{\mathrm{x}},\mathrm{BBN}}$ is not bimodal (compare with Figure˜4) perhaps due to a larger region explored by $\xi_{\nu_{e},\mathrm{BBN}}$, allowing $\xi_{\nu_{\mathrm{x}},\mathrm{BBN}}$ to remain closer to 0. When BBN is included, a bimodal appears for $\xi_{\nu_{\mathrm{x}},\,\mathrm{BBN}}$ for both the datasets. While for EMPRESS the bimodal is a signal of a preference for higher $N_{\mathrm{eff}}$, in the LBT case it might be due to not having properly explored the parameter region around $\xi_{\nu_{\mathrm{x}},\mathrm{BBN}}$ = 0, as the contour of LBT in the $Y_{\mathrm{He}}$-$\xi_{\nu_{\mathrm{x}},\mathrm{BBN}}$ lies at the 2$\sigma$ region. Particularly interesting, around recombination, we notice the appearance of a circumference in the plane of the two recombination nodes, $\xi_{\nu_{\mathrm{e}},\,\mathrm{rec}}$ and $\xi_{\nu_{\mathrm{x}},\,\mathrm{rec}}$, due to the level set of constant (and larger than 3.044, see Table˜4) $N_{\mathrm{eff}}$, expressed in terms of the two chemical potentials.

We find also in this case a longer tail in the posterior distribution of $H_{0}$, partly lost once we include the BBN data, a consequence of the tail on $N_{\mathrm{eff,\,BBN}}$. It is interesting to see that even when including LBT, which should recover the standard predictions, we still have a higher value of $H_{0}$, alleviating the $H_{0}$ tension.

Finally, regarding the Bayes factors, we see moderate and strong evidences againts $\Lambda$CDM. We conjecture that this is due to the large number of parameters and to the time variation of some of them, able to accommodate possible discrepancies in the data, as the different $\Omega_{\mathrm{m}}$ preferred by BAO and CMB data, for example. As a matter of fact, this could provide a better fit to some features that otherwise are not accommodated by other analyses. Notice also that we are not reaching a level of convergence strong enough to actually reach significant conclusions on the results of this case and due to degeneracies between the parameters, we could be dealing with statistical features more than physical ones.