How to relax the cosmological neutrino mass bound

TL;DR

This paper shows that the cosmological neutrino mass bounds are highly sensitive to the assumed neutrino momentum distribution, allowing for larger neutrino masses if their average momentum exceeds the thermal expectation.

Contribution

It demonstrates that non-standard neutrino momentum distributions can significantly relax cosmological neutrino mass bounds, highlighting a degeneracy with $N_{\mathrm{eff}}$ and $m_{\nu}$.

Findings

Neutrino distribution degenerates with $N_{\mathrm{eff}}$ and $m_{\nu}$ in observable effects.

Mass bounds can be relaxed by a factor of two with larger average neutrino momentum.

Current cosmological data cannot uniquely determine neutrino mass without assumptions on momentum distribution.

Abstract

We study the impact of non-standard momentum distributions of cosmic neutrinos on the anisotropy spectrum of the cosmic microwave background and the matter power spectrum of the large scale structure. We show that the neutrino distribution has almost no unique observable imprint, as it is almost entirely degenerate with the effective number of neutrino flavours, $N_{\mathrm{eff}}$, and the neutrino mass, $m_{\nu}$. Performing a Markov chain Monte Carlo analysis with current cosmological data, we demonstrate that the neutrino mass bound heavily depends on the assumed momentum distribution of relic neutrinos. The message of this work is simple and has to our knowledge not been pointed out clearly before: Cosmology allows that neutrinos have larger masses if their average momentum is larger than that of a perfectly thermal distribution. Here we provide an example in which the mass limits are relaxed by a factor of two.

Tables

Table 1. Table 1 : 95% CL upper bounds on the sum of neutrino masses for different choices of the multipoles { F 1 , F 2 } subscript 𝐹 1 subscript 𝐹 2 \{F_{1},F_{2}\} ( 2.8 ) and the normalization N 𝑁 N fixed to ( 3.2 ), in comparison to the standard FD distribution. All upper limits are based on Planck 2015 likelihoods [ 8 ] .

	TT+lowP (95 % CL)	TT+lowP+BAO (95 % CL)
FD	eV	eV
$F_1=0.92,F_2=0$	eV	eV
$F_1=0,F_2=1.06$	eV	eV
$F_1=0.92,F_2=1.06$	eV	eV

Equations

$${\dot{\Psi }}_{\nu ,0}$$

$${\dot{\Psi }}_{\nu ,2}$$

$$\delta {\rho }_{\nu }$$

$$({\bar{\rho }}_{\nu }+{\bar{P}}_{\nu }){\sigma }_{\nu }$$

$$f_{\nu }(x)=\frac{1}{e^x+1}.$$

$${\int }_0^{\infty }\mathrm{d}x\frac{1}{e^x+1}{p}_n(x)p_m(x)={\delta }_{nm}.$$

$$f_{\nu }(x)=\frac{1}{e^x+1}\sum _{n=0}^{\infty }{c}_n{p}_n(x),$$

$$c_n={\int }_0^{\infty }\mathrm{d}x{f}_{\nu }(x)p_n(x).$$

$$f_{\nu}(x)=N\cdot\frac{1}{e^{x}+1}\bigg{(}p_{0}(x)+F_{1}p_{1}(x)+F_{2}p_{2}(x)\bigg{)},$$

$$0\le F_1<0.92,0\le F_2<\mathrm{1.06.}$$

$$f_{\nu }^{\mu }(x)=\frac{1}{e^{x+\mu }+1},$$

$$N=c_0=\mathrm{0.665848.}$$

$$f_{\nu }^y(x)=\frac{1}{e^x+1}\left[1+y\frac{x{e}^x}{e^x+1}\left(x\frac{e^x-1}{e^x+1}-4\right)\right].$$

$$N=c_0=0.499533,F_1=\frac{c_1}{N}=7.921424\cdot 10^{-8},F_2=\frac{c_2}{N}=\mathrm{0.775381.}$$

$$N\cdot\int\mathrm{d}x\frac{x^{3}}{e^{x}+1}\bigg{(}p_{0}(x)+F_{1}p_{1}(x)+F_{2}p_{2}(x)\bigg{)}\overset{!}{=}\int\mathrm{d}x\frac{x^{3}}{e^{x}+1},$$

$$N=\frac{5.682197}{6.825012+18.239415F_1+17.420046F_2}.$$

$$\delta C_{\ell }=\frac{C_{\ell }-C_{\ell }^{FD}}{C_{\ell }^{FD}}$$

$$\int_{0}^{\infty}\mathrm{d}x\,x^{2}\sqrt{x^{2}+\frac{m_{\nu}^{2}}{T_{\nu 0}^{2}}}\,\left(\frac{1}{e^{x}+1}\right)\overset{!}{=}N\int_{0}^{\infty}\mathrm{d}x\,x^{2}\sqrt{x^{2}+\frac{{m_{\nu}^{*}}^{2}}{T_{\nu 0}^{2}}}\,\frac{1}{e^{x}+1}\bigg{(}p_{0}(x)+F_{1}p_{1}(x)+F_{2}p_{2}(x)\bigg{)},$$

$$F_1$$

$$p_0(x)$$

$$p_1(x)$$

$$p_2(x)$$

$$p_3(x)$$

$$p_4(x)$$

$$p_5(x)$$

$$p_0(x)+F_1{p}_1(x)+F_2{p}_2(x),$$

Full text

1 Introduction

The neutrino sector is one of the most obscure sectors within the standard model of particle physics. The nature of neutrinos, the neutrino mass production mechanism and even the values of the neutrino masses are still open puzzles and matter to an active field of research.

Since the observation of neutrino oscillations revealed that neutrinos are massive, the determination of the neutrino masses has become the holy grail of neutrino physics. While the neutrino mass splittings have already been measured [1, 2], the neutrino mass ordering as well as the mass of the lightest mass state are still undetermined. The strongest bound on the neutrino mass from laboratory searches currently comes from the Troitzk and Mainz experiments, namely $m_{\nu_{e}}<2.05$ eV (95% CL) [3] and $m_{\nu_{e}}<2.3$ eV (95% CL) [4], where $m_{\nu_{e}}=(\sum_{i}|U_{ei}|^{2}m_{i}^{2})^{1/2}$ and $U_{ei}$ are the entries of the Pontecorvo–Maki–Nakagawa–Sakata matrix. The tritium decay experiment KATRIN aims to reach a sensitivity of $m_{\nu_{e}}<$0.2 eV at 90% CL [5].

The by far strongest constraint on the sum of the neutrino masses today however comes from cosmological observations. Measurements of the temperature and polarization anisotropy spectrum of the cosmic microwave background (CMB) with the Planck satellite in combination with other cosmological measurements like baryonic acoustic oscillations (BAO) show that $\sum m_{\nu}<0.12$ eV (95% CL) [6]. This tight constraint has also been derived in the analysis of [7] including galaxy power spectra but with Planck 2015 data [8]. Future data analyses promise even tighter bounds and will be realized in the next decade [9, 10, 11, 12, 13, 14, 15]. The cosmological mass bounds relax to some extent when marginalized over different extensions of the standard $\Lambda$CDM model [17] but were found to be relatively robust [8, 16].

While most of the $\Lambda$CDM extensions studied in [8, 11, 16, 9, 14, 12, 17] are addressing non-neutrino properties, we show how the cosmological neutrino mass bound can be lowered by directly targeting an assumption that we usually make about cosmological neutrinos: We abandon the assumption of a Fermi-Dirac distribution for neutrinos and allow non-thermal features in the neutrino phase-space distribution function. The idea behind is as simple as follows: By allowing relic neutrinos to have a higher average-momentum they can have much larger masses.

We also address the question which impact the neutrino distribution function has in the first place on the CMB and the large scale structure (LSS) of the Universe. We show that there is almost no unique feature of the neutrino distribution function in the CMB and LSS since non-thermal distortions can almost be entirely mimicked by changes in the values of the effective number of neutrino flavours $N_{\text{eff}}$ and $\sum m_{\nu}$.

While the degeneracy between the neutrino distribution function and $N_{\text{eff}}$ is widely known, its intrinsic degeneracy with $\sum m_{\nu}$ has been much less communicated. The latter has however already been pointed out in [18]. But at the time of this work cosmological data were still much less precise than today. In particular at the time of [18, 19, 20], there was still a large degeneracy between $\sum m_{\nu}$ and $N_{\mathrm{eff}}$ which also allowed to significantly lower the cosmological neutrino mass bound. Due to today’s precise CMB data of the Planck satellite [8] the degeneracy between $N_{\mathrm{eff}}$ and $\sum m_{\nu}$ is mainly eliminated such that a non-standard $N_{\mathrm{eff}}$ alone does not significantly weaken the neutrino mass bound anymore. The impact of a non-thermal neutrino distribution function on the neutrino mass bound however has to our knowledge not been studied after the first release of Planck data in 2013 [21].

This work is organized as follows: In section 2 we present a parameterization for deviations from equilibrium of the neutrino distribution function. It describes out-of-equilibrium distributions that are not too far from the thermal case. As examples, we apply our parameterization to non-standard distributions inspired by $\mu$- and $y$-type distortions of the cosmic photon spectrum [22]. We present the impact of our parameterization on the CMB temperature anisotropy spectrum and the LSS matter power spectrum in section 3. By normalizing the distribution function and adopting an effective neutrino mass, we show that our parameterization allows to weaken the neutrino mass bound by more than 100%. We explicitly demonstrate this by performing a Markov chain Monte Carlo (MCMC) analysis in section 4 and finally conclude in section 5.

2 Relic neutrino distribution

For massless neutrinos cosmological perturbation theory is entirely independent of the neutrino distribution function. All relevant quantities related to neutrinos – e.g. the energy contrast or velocity divergence – are already integrated over momentum. However, for massive neutrinos the Boltzmann hierarchy is momentum-dependent [23]. The equations for the zeroth and second multipole moments ($\Psi_{\nu,\ell}$) directly depend on the neutrino background-distribution function $f_{\nu}$, i.e. (following the notation of [23])

[TABLE]

Here, $q$ denotes comoving momentum and $\epsilon$ comoving neutrino energy, $k$ wave number, $h$ and $\eta$ are metric perturbations in the synchronous gauge, and overdots denote derivatives with respect to conformal time. The contributions to the energy-momentum tensor are derived by solving (2.1) together with the $\ell=1$ and $\ell>2$ equations on a momentum grid, and integrating over momentum, e.g.

[TABLE]

for the energy density perturbation $\delta\rho_{\nu}$ and the anisotropic stress $\sigma_{\nu}$. $a$ denotes the scale factor. Hence, the neutrino distribution function $f_{\nu}$ has a direct impact on the evolution of massive neutrino perturbations which in turn have an impact on the photon and matter perturbations. Cosmological observables like the CMB anisotropy spectrum and the LSS matter power spectrum therefore depend on the neutrino distribution as well.

Let us discuss how to parameterize non-standard neutrino distributions in a relatively model-independent way in the the following section.

2.1 Parameterization of relic neutrino distribution

According to the standard model of particle physics, neutrinos decouple from the cosmic plasma at around $T\sim$ 1 MeV. Since neutrinos are effectively massless at decoupling, they are usually assumed to follow an ultra-relativistic Fermi-Dirac (FD) distribution, which is even preserved as the neutrino momentum distribution at later times when neutrinos become non-relativistic,

[TABLE]

Here, we introduced the dimensionless comoving momentum $x=q/T_{\nu,0}$, where $T_{\nu,0}$ denotes the neutrino temperature today. Note that there are small distortions to the FD distribution (2.4) due to $e^{+}-e^{-}$ annihilation which are however usually expressed in terms of an increased number of neutrino flavours, i.e. $N_{\text{eff}}=3.046$ [24, 25].

There are many exotic scenarios which could possibly alter the standard distribution function (2.4) of neutrinos. Examples for this are: the decay or annihilation of some heavier particle into neutrinos (e.g. [26, 27, 18, 28]), non-standard neutrino interactions like the freeze-in of a massless Majoron (e.g. [29, 30]) or decaying neutrinos (e.g. [31, 32])111Non-standard interactions scenarios in general would however also change the free-streaming behavior of neutrinos [33] which would induce further modifications to the CMB and LSS that are not discussed here., or alternatively oscillations to sterile neutrino states [34, 35, 36, 37, 38]. In order to keep the discussion as general as possible we do not want to focus on any specific scenario but aim to be as model-independent as possible. We therefore choose a parameterization very similar to [39] and expand the distribution function in polynomials which are orthonormal with respect to the FD distribution, i.e.

[TABLE]

This choice for the orthonormality relation is different from [18], where the polynomials were chosen to be orthonormal to $x^{2}/(e^{x}+1)$. Since one of the aims of this work is to study the impact of the distribution function itself (i.e. of $f_{\nu}(x)$ and not of $x^{2}f_{\nu}(x)$), we find our choice more appropriate for this purpose. The first six polynomials can be found in the appendix A. The distribution function can then be written as

[TABLE]

where the coefficients $c_{n}$ are given by

[TABLE]

Note that equation (2.6) is entirely model-independent since any function can be decomposed into the orthogonal polynomials (2.5). In general however, a distribution function very different from the thermal FD distribution necessitates the inclusion of a high number of coefficients $c_{n}$ (2.7). Nevertheless, by restricting our discussion to only a few moments in (2.6), we can already describe a relatively large variety of distribution functions. Let us therefore consider only polynomials up to order $n=2$ in the following.

A disadvantage of the parameterization (2.6) is that for arbitrary choices of $c_{0}$, $c_{1}$ and $c_{2}$ the distribution can get negative values which is of course an unphysical feature. It is however easy to show that this unphysical behavior can be avoided by restricting the multipoles to $c_{1}\leq 0.92c_{0}$ and $c_{2}\leq 1.06c_{0}$ (see appendix A for a detailed discussion). In order to better control $c_{1}$ and $c_{2}$ it turns out to be more convenient to factorize $c_{0}$ in (2.6) and interpret it as an overall normalization $N$. Without loss of generality, the expansion (2.6) up to $n=2$ can therefore be rewritten as

[TABLE]

where $N\equiv c_{0}$ and $F_{1/2}\equiv(c_{1/2}/N)$. The conditions to have a positive distribution function therefore translate into

[TABLE]

2.2 Examples

As an example, we apply our expansion (2.8) to two types of familiar non-thermal features in distribution functions, namely $\mu$- and $y$-type distortions. Both of these distortions are known from the CMB spectrum, where they can appear due to inefficient ($y$-type) or efficient ($\mu$-type) Compton scattering with electrons [22]. $\mu$-type distortions could furthermore be realized by neutrino chemical potentials, which would appear if the lepton asymmetry of our Universe was non-negligible, see e.g. [40, 41]. In that case, they should however not be interpreted as distortions since a chemical potential is a thermal (but non-standard) feature of the neutrino distribution function. For $y$-type distortions there is however no obvious realization within the neutrino sector. In analogy to the Sunyaev-Zeldovich effect [22], such a distortion could in principle be produced by neutrinos scattering on some high-energetic particles after neutrino decoupling. However, without claiming any specific production mechanism of such a distortion, we consider it – only for the purpose of testing how well the parameterization (2.8) can mimic non-standard distribution functions.

Let us start with considering a $\mu$-type distortion,

[TABLE]

where we exemplary chose $\mu=0.3$ here. In thermodynamic equilibrium $\mu$ denotes the ratio of chemical potential and temperature and after neutrino decoupling it becomes the so-called neutrino degeneracy parameter.

The normalization $N$ in (2.8) can be calculated according to (2.7) and is given by

[TABLE]

The $F_{1}$ and $F_{2}$ parameters are much smaller than 1, which shows that they are very subdominant in comparison to the overall normalization $N$. Indeed, if we neglect $F_{1}$ and $F_{2}$ (i.e. we cut the expansion at $n>0$), we can see in the left panel of figure 1 that the distribution function resembles the exact expression in (2.10) quite well.

$y$-type distortions are described according to

[TABLE]

Note that compared to $y$-type distortions of CMB photons there are some sign switches in (2.12) due to the difference between bosons and fermions in the Kompaneets equation.

Exemplary we study the case of $y=0.2$ and find for the expansion parameters in (2.8)

[TABLE]

The right panel of figure 1 shows the exact expression (2.12) with $y=0.2$ and the expansion (2.8) with the parameter values (2.13) in comparison. Even though the expansion does not accurately reproduce the exact expression near the local minimum, the overall consistency is reasonably good.

We conclude that our expansion can mimic $\mu$- and $y$-type distortions. However, as we include only quadratic polynomials in equation (2.8) our parameterization is of course restricted to certain functions and we cannot mimic every thinkable distribution function by equation (2.8). Given a specific model for non-standard neutrino distributions it may be more useful to parameterize the distribution function in a different way.

This being said, we also want to point out another conclusion that we can draw from this section: Since the normalization $N$ has a very similar effect as the effective number of neutrino flavours $N_{\mathrm{eff}}$, we conclude that chemical potentials are degenerate with $N_{\mathrm{eff}}$. While this is a known fact for massless neutrinos (where the distribution function itself does not matter, as explained above), this has not been a priori clear for the case of massive neutrinos. Therefore, even with future experiments like CMB-S4 [42] or Euclid [43, 44] it will be impossible to disentangle the effect of neutrino chemical potentials and $N_{\mathrm{eff}}$. Note that this conclusion however only holds for the muon and tau neutrino, since the electron neutrino chemical potential has an additional impact on big bang nucleosynthesis (BBN) that changes the primordial helium abundance and has a much larger impact on the CMB anisotropy spectrum [41].

3 Impact of neutrino spectral distortions on the CMB and LSS

Let us now study in detail the impact of the parameterization (2.8) of neutrino spectral distortions on the CMB temperature anisotropy spectrum and the LSS matter power spectrum. In order to do so we should make a short detour to BBN.

The impact of neutrinos on BBN is two-fold: The electron-neutrino distribution function directly impacts BBN as it determines the neutron-to-proton ratio at the onset of BBN at $T=\mathcal{O}$(0.1 MeV). The impact of the muon- and tau-neutrinos is more indirect as they only impact BBN through the expansion rate, which depends on the energy density of the ultra-relativistic neutrinos and hence can be parametrized by $N_{\text{eff}}$. Hence, in principle a non-standard neutrino distribution function changes the BBN prediction for the abundance of primordial helium, which in turn has an impact on the CMB anisotropy spectrum. However, in most thinkable scenarios the neutrino distribution function would probably still be FD at BBN: Before neutrino decoupling at $T=\mathcal{O}$(1 MeV) any non-thermal feature in the neutrino distribution function would be washed out quickly by weak interactions 222Note that this does not hold for chemical potentials – discussed in the previous section – which are a thermal (but non-standard) feature in the neutrino distribution function.. The window for non-standard physics to alter the neutrino distribution before the onset of BBN is therefore restricted to the cosmic temperature range $T=\mathcal{O}$(0.1 MeV) to $T=\mathcal{O}$(5 MeV). For a detailed analysis of such a scenario see [26, 27]. However, in the general framework of this work it seems to be reasonable to restrict our discussion to scenarios where the neutrino distribution only becomes non-standard after the onset of BBN. Some scenarios could nevertheless still modify BBN through an enhanced $N_{\text{eff}}$, e.g. the decay of a particle of keV mass. Note that a particle of $>\mathcal{O}$(MeV) mass would in contrast only enhance $N_{\mathrm{eff}}$ after BBN, as it would already be non-relativistic and therefore not enhance the relativistic energy density at BBN. As we discuss in the next section, we focus on scenarios that do not alter $N_{\mathrm{eff}}$ but only the way how neutrinos are distributed.

3.1 Degeneracy of $N_{\text{eff}}$ and spectral distortions — normalization of the distribution function

It is already well-known that the neutrino distribution function is degenerate with $N_{\mathrm{eff}}$, see e.g. [18]. In order to study the required sensitivity of current or future experiments to disentangle the effect of $N_{\mathrm{eff}}$ and the distribution function, we are interested in finding out if there is a unique feature of the neutrino distribution function in the CMB anisotropy spectrum and the LSS matter power spectrum. In order to switch off the effect due to $N_{\mathrm{eff}}$ it is useful to choose the normalization $N$ (2.8) such that $N_{\text{eff}}$ always gives the same value. This can be realized if

[TABLE]

which is fulfilled by the following choice of $N$

[TABLE]

Such a scenario (i.e. a non-thermal distribution function but an unaltered $N_{\mathrm{eff}}$) could be realized for example by the freeze-in of a massless Majoron [29, 30], decaying neutrinos [31, 32] or oscillations to sterile neutrinos [34, 35, 37].

In figure 2 we show the normalized (3.2) neutrino distribution function (2.8) as well as $x^{3}f_{\nu}(x)$ for different choices of $F_{1}$ and $F_{2}$. All shown cases of spectral distortions give rise to an increase of the mean neutrino momentum as compared to the FD spectrum. We can see that the dotted line with $\{F_{1}=0,F_{2}=1.06\}$ is the most efficient parameter combination in shifting the average momentum.

3.2 Implementation and cosmological observables

We implemented the parameterization (2.8) with the normalization given in (3.2) into the linear Einstein-Boltzmann solver class [45]. Note that for parameter combinations with $F_{1}\neq 0$ and $F_{2}\neq 0$ at the same time the precision of the numerical momentum-integration method has to be slightly increased in order to obtain precise results. Here we focus our study on mass degenerate neutrinos, i.e. $m_{1}=m_{2}=m_{3}$, as we are only interested in the upper limit on the sum of all neutrino masses. The base model is a spatially flat, $\Lambda$ cold dark matter model with scale-free primordial curvature fluctuations and without primordial gravitational waves, i.e. todays vanilla model.

Figure 3 shows the relative percent difference in the angular power spectrum of CMB temperature anisotropy ($C_{\ell}$) and the LSS matter power spectrum ($P(k)$) between a non-standard neutrino distribution function (2.8) and the standard FD distribution, i.e.

[TABLE]

and equivalently for the matter power spectrum $P(k)$. The dark blue curve ($\sum m_{\nu}=0$ eV) is only plotted as a numerical test to confirm that for massless neutrinos cosmological observations are indeed independent of the distribution function.

As expected, the impact of a non-thermal distribution function is increasing for increasing neutrino masses (in fact approximately proportional to the sum of neutrino masses). Remarkably the relative difference in the LSS matter power spectrum is up to 20% for the parameter combinations under consideration, whereas it is about a factor of 3 smaller – but still significant – in the CMB anisotropy spectrum. Besides a modulation of the angular power spectrum, both in the CMB and in the LSS we can observe an enhancement of structure at small length scales. This shows that the neutrino distribution function has indeed an impact on the the CMB and LSS which is not degenerate with $N_{\mathrm{eff}}$. The remaining question is now whether the effects shown in figure 3 are unique features of the neutrino distribution function or if they can be mimicked by other cosmological parameters.

3.3 Determining an effective neutrino mass

Of course, the normalization $N$ (3.2) that we chose in the last section fixes the neutrino energy density at early times when neutrinos are relativistic. There is no way to find a normalization such that the energy density of the non-thermal distribution (2.8) is always the same as in the standard FD case. However, the energy density at late times is not only determined by the distribution function but also by the neutrino mass. The CMB anisotropy spectrum of figure 3 also shows that the impact of a given neutrino mass can be mimicked to some extent by a different neutrino mass and another corresponding distribution function: For example, the red solid line ($\sum m_{\nu}=0.6$ eV, $F_{1}=0.92$, $F_{2}=0$) looks very similar to the green dotted line ($\sum m_{\nu}=0.3$ eV, $F_{1}=0.92$, $F_{2}=1.06$). This motivates the question if we can eliminate the signal of a non-thermal neutrino distribution in figure 3 by adapting the neutrino mass accordingly. In other words, we want to determine the neutrino mass of the non-thermal distribution $m^{*}_{\nu}$, such that today’s $\rho_{\nu}(m^{*}_{\nu})=\rho^{FD}_{\nu}(m_{\nu})$. Then the following equation is fulfilled

[TABLE]

where $N$ is given in equation (3.2).

By approximately solving equation (3.4) we find the following relations for the parameter combinations used in figure 3:

[TABLE]

Figure 4 shows how the signal presented in figure 3 changes when adapting the neutrino masses according to (3.5). After adapting the neutrino masses the remaining effect of the spectral distortions drops below the 1% level in both, temperature anisotropy spectrum and LSS matter power spectrum. Note that the numerical precision of our calculations are at the level of $0.02\%$ for the CMB anisotropy spectrum and about $0.2\%$ for the LSS matter power spectrum (dark blue curves). In order to track down how much of the signal in figure 4 is due to the real unique feature of the neutrino distribution function and how much is simply numerical noise, we would have to tune the numerical precision of class. However, the remaining signal in fig. 4 is anyhow below the precision of current experiments and its detection therefore currently unfeasible.

We have just shown that the neutrino distribution function is not only degenerate with $N_{\mathrm{eff}}$ but also with the neutrino mass. This has interesting consequences on the robustness of cosmological neutrino mass bounds: By abandoning the assumption that relic neutrinos follow a thermal FD distribution precisely, we can weaken cosmological neutrino mass bounds. As we have shown, cosmological observations almost exclusively care about the time when most of the neutrinos become non-relativistic. This time is not only determined by the neutrino mass but also by the neutrino distribution function. Even when keeping the neutrino energy density at early times fixed (i.e. when $N_{\mathrm{eff}}$ is fixed), we can still enhance the neutrino average momentum such that neutrinos are allowed to be much heavier than in the standard framework. For the parameterization chosen in this work (2.8) we expect the neutrino mass bounds to become weaker by more than $90\%$ (for $\{F_{1}=0,F_{2}=1.06\}$, as shown in figure 4). We however want to strongly emphasize here that the amount by which the neutrino mass bound is changed clearly depends on the specific model under consideration. Some models with a non-thermal neutrino distributions may lead to a weaker effect, while other models may even relax the neutrino mass bound by much more than the $90\%$ found here.

4 MCMC analysis

In order to demonstrate on real data that we can lower the neutrino mass bound by allowing a non-thermal neutrino distribution, we perform an MCMC analysis with the MCMC engine Montepython [46, 47]. We study the two following data sets representing conservative and more stringent cosmological mass bounds:

TT+lowP

Temperature anisotropy spectrum and low-$\ell$ polarisation from the Planck 2015 data [8].

TT+lowP+BAO

Same as above plus measurements of the baryon acoustic oscillation peak scale as measured by 6DF [48], BOSS LOWZ and CMASS [49, 50], and the SDSS Main Galaxy Sample [51].

As in the previous section, we study the case of mass-degenerate neutrinos since the neutrino mass-splittings are still close enough to be negligible with current data. We apply flat priors on all cosmological parameters and adopt the Gelman-Rubin convergence criterion $R<0.01$.

We study different distribution functions, where we choose the same combinations of the expansion parameters $\{F_{1},F_{2}\}$ (2.8) as in figure 4 and equation (3.5). The overall normalization $N$ is thereby again fixed according to equation (3.2). This also guarantees that we do not introduce any inconsistency with the CMB constraints on $N_{\text{eff}}$, which is fixed to $N_{\text{eff}}=3.046$ in this analysis. Note that a global fit, where $F_{1}$ and $F_{2}$ and $m_{\nu}$ are treated as free parameters, does not provide any further insight: Since at $m_{\nu}=0$ cosmic perturbation theory is blind to the relic neutrino distribution, $F_{1}$ and $F_{2}$ are entirely unconstrained in such an analysis. This also implies that future cosmological data will only have constraining power on the relic distribution function if they turn out to exclude $\sum m_{\nu}=0$ or if the neutrino masses eventually will be measured by laboratory experiments.

As already noted in section 3, the precision of the numerical integration method for massive neutrinos in class needs to be enhanced in order to properly perform momentum integrations. As a comparison we also performed an MCMC run assuming the standard FD distribution. Our results for the mass constraints in this run agree well with the ones by the Planck collaboration [8], which are however weaker than the ones quoted in the introduction, since the Planck 2018 data [6] are not publicly available at the time of finishing this work.

Figure 5 shows the posterior distribution of the neutrino mass sum $\sum m_{\nu}$ for the different parameter combinations. We clearly see how the posterior stretches, depending on the assumptions on the neutrino momentum distribution. The 95% CL upper bounds for the different parameter combinations and different data sets are reported in table 1. The mass bounds we find in our analysis confirm the prediction of equation (3.5), i.e. depending on the choices of the parameters $\{F_{1},F_{2}\}$ we can lower the bound by more than 90%, in fact by even more than a factor of two in the case of the more stringent bounds (0.37 eV instead of 0.18 eV).

5 Conclusions

In this work, we studied the impact of a non-thermal relic neutrino distribution on cosmological observables, i.e. the CMB temperature anisotropy spectrum and the LSS matter power spectrum. Any process that increases the mean momentum of cosmic relic neutrinos could give rise to an increased cosmological upper bound on the sum of the neutrino masses. Two types of processes can be imagined. Either a process that modifies the bulk of cosmic relic neutrinos, or a rare effect that affects only a small number of cosmic neutrinos but injects a lot of energy, e.g. a post-BBN decay of a heavy particle into high-energy neutrinos. Our formalism is well suited to study the former type of effects, but the latter one requires much more care with the numerical accuracy of the Boltzmann codes to account for the hierarchy in neutrino energy scales. Here we have focused on the former type of effects.

In order to parameterize non-thermal distributions in a model-independent way, we expanded the distribution function in orthonormal polynomials where we cut the expansion at order $n\geq 3$ (2.8) for simplicity. This introduces three parameters, $N$, $F_{1}$ and $F_{2}$. We demonstrated that our parameterization can mimic $\mu$-distortions and $y$-inspired distortions very well. In order to figure out if the distribution function has a unique imprint in the CMB and the LSS we chose the normalization $N$ such that it always reproduces the standard value of $N_{\mathrm{eff}}$. A realization of such a scenario would be any (non-standard) physics that redistribute neutrinos without injecting a too large amount of energy into the neutrino sector. Examples of this are e.g. non-standard neutrino interactions like the freeze-in of a massless Majoron (e.g. [29, 30]), decaying neutrinos (e.g. [31, 32]), or oscillations of active to sterile neutrinos [34, 35, 36, 37, 38]. After normalizing the distribution function, the remaining effect is up to 6% in the CMB temperature anisotropy spectrum, and up to 20% in the LSS matter power spectrum, depending on the choices of $\{F_{1},F_{2}\}$ and $\sum m_{\nu}$. For different choices of the parameters $\{F_{1},F_{2}\}$ we furthermore calculated an effective neutrino mass $m_{\nu}^{*}$ that allows to also have the same energy density at late times as the standard FD distribution with mass $m_{\nu}$, see equations (3.4) and (3.5). When adjusting the neutrino mass of the non-thermal distribution to the effective mass $m_{\nu}^{*}$ the signal in the CMB and LSS shrinks down below the percent level. This means that the class of spectral distortions studied in this work has almost no unique observable imprint, since it is almost entirely degenerate with $N_{\mathrm{eff}}$ and $\sum m_{\nu}$. While the degeneracy to $N_{\mathrm{eff}}$ has already been well-known, the intrinsic degeneracy to $\sum m_{\nu}$ has been much less communicated and — in our opinion — overseen in discussions about the robustness of cosmological neutrino mass bounds. By performing an MCMC analysis we have demonstrated that the cosmological neutrino mass bound can be relaxed by more than 90% for the parameterization studied in this work. Note however that our parameterization (2.8) is constructed in such a way that it only describes distribution functions that are still relatively close to the standard FD distribution. Adding higher polynomials to the parameterization could further relax the mass bounds and assuming a drastically different distribution function could in principle relax the mass bounds in a drastic way. In general, given a specific model it has to be tested from case to case if our parameterization (2.8) is sufficient in mimicking the neutrino distribution and by which amount the neutrino mass bound would be relaxed.

To sum up, we showed how much the cosmological neutrino mass bound relies on our assumption about the relic neutrino distribution function. This stands in contrast to claims that the mass bound is entirely robust. Indeed, assuming that neutrinos follow a thermal FD distribution is — on theoretical grounds — reasonable. But from an observationally agnostic point of of view, we have never measured the relic neutrino momentum distribution (and will not in the near future), so we should keep in mind that the cosmological mass bound and the constraints on $N_{\text{eff}}$ strongly depend on this assumption. Direct mass detection experiments like KATRIN [5] are therefore still of great importance, even though their sensitivity to the neutrino mass is already beaten by cosmological observations under the assumption of thermally distributed relic neutrinos. If KATRIN measures a neutrino mass in contradiction to the cosmological neutrino mass bound, this would have exciting consequences for neutrino physics and the most promising cosmological model would probably be provided by scenarios that alter the momentum distribution of cosmic relic neutrinos. Conversely, an improved upper limit on the electron neutrino mass from KATRIN would increase our confidence in the thermal history of the early Universe.

Acknowledgements

We acknowledge financial support by MEC and FEDER (EC) Grants No. SEV-2014-0398, FIS2015-2245-EXP, FPA2014-54459 and the Generalitat Valenciana under Grant No. PROME-TEOII/2013/017 (IMO and GB), by the Deutsche Forschungemeinschaft (DFG, German Research Foundation) Research Training Group 1620 ‘Models of Gravity’ (SK and DJS) and Project number 315477589 TRR 211 (IMO and DJS) and by EU Networks FP10 ITN ELUSIVES (H2020-MSCA-ITN-2015-674896) and INVISIBLES-PLUS (H2020-MSCA- RISE-2015-690575) (JS and GB), and MINECO grant FPA2016-76005- C2-1-P and Maria de Maetzu grant MDM-2014-0367 of ICCUB (JS). Numerical calculations have been performed using the computing resources from the Tirant cluster of University of Valencia.

Appendix A Orthogonal Polynomials

The first six polynomials fulfilling the orthonormality relation (2.5) are given by

[TABLE]

An illustration of the polynomials (A.1)-(A.6) can be found in figure 6. In the following, we however only focus on the first three polynomials.

As discussed in section 2.1, the expansion in those polynomials (2.8) has the disadvantage of evolving negative values for the distribution function for certain choices of parameter combinations $\{F_{1},F_{2}\}$. This unphysical feature is caused by the orthogonal polynomials themselves having negative values in certain ranges of comoving momentum. In particular, the term in brackets in equation (2.8),

[TABLE]

can become negative.

Let us first study the scenario of $\{F_{1}\neq 0,F_{2}=0\}$. The $p_{1}(x)$ polynomial (A.2) has its global minimum at $x=0$. That implies that as long as we have $F_{1}<-p_{0}(0)/p_{1}(0)\approx 0.92$, equation (A.7) and the distribution function (2.8) are positive for all $x$. The lower limit (i.e. $F_{1}>0$) is derived by demanding the correct asymptotic behaviour of $f_{\nu}(x)$ (2.8) for $x\rightarrow\infty$.

For $\{F_{1}=0,F_{2}\neq 0\}$ the same argument gives $F_{2}<-p_{0}(x_{0})/p_{2}(x_{0})\approx 1.06$, where $x_{0}=2.141462$ is the global minimum of $p_{2}(x)$. For mixed states, i.e. $\{F_{1}\neq 0,F_{2}\neq 0\}$, the condition of a positive distribution function for $F_{2}$ depends on $F_{1}$ (and the one for $F_{1}$ on $F_{2}$). On general grounds (particularly for an MCMC analysis), such conditions are less convenient to handle. For the sake of simplicity, we therefore hold on the conditions (2.9), even though we thereby omit some of the possible parameter space of $\{F_{1},F_{2}\}$. By restricting our considerations to the parameter ranges (2.9), we however still cover a large variety of distribution functions.