Cosmic microwave background constraints on secret interactions among sterile neutrinos

TL;DR

This paper uses Planck CMB and BAO data to constrain secret interactions among sterile neutrinos, finding significant tension with cosmological observations and setting limits on neutrino mass and interaction strength.

Contribution

It provides the first comprehensive cosmological constraints on secret sterile neutrino interactions, including their impact on perturbations and mass bounds.

Findings

Sterile neutrino mass bound: m_s < 0.82 eV (Planck alone)

Interaction strength limit: G_X < 2.8 × 10^{10} G_F

Disfavors models with M_X ~ 0.1 MeV and g_X ~ 0.1

Abstract

Secret contact interactions among eV sterile neutrinos, mediated by a massive gauge boson $X$ (with $M_X \ll M_W$), and characterized by a gauge coupling $g_X$, have been proposed as a mean to reconcile cosmological observations and short-baseline laboratory anomalies. We constrain this scenario using the latest Planck data on Cosmic Microwave Background anisotropies, and measurements of baryon acoustic oscillations (BAO). We consistently include the effect of secret interactions on cosmological perturbations, namely the increased density and pressure fluctuations in the neutrino fluid, and still find a severe tension between the secret interaction framework and cosmology. In fact, taking into account neutrino scattering via secret interactions, we derive our own mass bound on sterile neutrinos and find (at 95% CL) $m_s < 0.82$ eV or $m_s < 0.29$ eV from Planck alone or in combination with BAO, respectively. These limits confirm the discrepancy with the laboratory anomalies. Moreover, we constrain, in the limit of contact interaction, the effective strength $G_X$ to be $ < 2.8 (2.0) \times 10^{10}\,G_F$ from Planck (Planck+BAO). This result, together with the mass bound, strongly disfavours the region with $M_X \sim 0.1$ MeV and relatively large coupling $g_X\sim 10^{-1}$, previously indicated as a possible solution to the small scale dark matter problem.

Tables

Table 1. Table 1 : Description of the models considered in this work.

	Description
$\Lambda \mathrm{CDM}$	Standard six-parameter $\Lambda \mathrm{CDM}$, $N_{\mathrm{eff}}=3.046$.
S$\Lambda \mathrm{CDM}$_$\mathrm{GX0}$	Sterile neutrino extension, $N_{\mathrm{eff}}=2.7$, $m_s$ free, “small” $G_X$ ($\sim {10}^8{G}_F$).
S$\Lambda \mathrm{CDM}$	Sterile neutrino extension, $N_{\mathrm{eff}}=2.7$, $m_s$ and $G_X$ free.
S$\Lambda \mathrm{CDM}$_Narrow	Sterile neutrino extension, $N_{\mathrm{eff}}=2.7$, $G_X$ free, $m_s=1.27\pm 0.03\mathrm{eV}$ (gaussian prior).
S$\Lambda \mathrm{CDM}$_Broad	Sterile neutrino extension, $N_{\mathrm{eff}}=2.7$, $G_X$ free, $0.93\mathrm{eV}\le m_s\le 1.43\mathrm{eV}$ (flat prior).

Table 2. Table 2 : Parameter constraints for the models under consideration, from the PlanckTT+lowP dataset. We either quote constraints in the form “mean ± plus-or-minus \pm 68% uncertainty”, or as 95% credible intervals (when not indicated, the lower limit should be understood to be zero). Units of m s subscript 𝑚 𝑠 m_{s} and H 0 subscript 𝐻 0 H_{0} are eV and km s -1 Mpc -1 , respectively.

Parameter	$\Lambda \mathrm{CDM}$	S$\Lambda \mathrm{CDM}$_GX0	S$\Lambda \mathrm{CDM}$	S$\Lambda \mathrm{CDM}$_Broad	S$\Lambda \mathrm{CDM}$_Narrow
${\Omega }_b{h}^2$	$0.02222\pm 0.00023$	$0.02177\pm 0.00024$	$0.02172\pm 0.00025$	$0.02167\pm 0.00025$	${0.02166}_{-0.00024}^{+0.00024}$
${\Omega }_c{h}^2$	$0.1197\pm 0.0021$	$0.1167\pm 0.0022$	$0.1171\pm 0.0023$	$0.1165\pm 0.0022$	$0.1160\pm 0.0021$
$100{\theta }_{MC}$	$1.04085\pm 0.00047$	$1.04103\pm 0.00050$	${1.04323}_{-0.00073}^{+0.00091}$	$1.04319\pm 0.00074$	${1.04307}_{-0.00077}^{+0.0010}$
$\tau$	$0.078\pm 0.019$	$0.070\pm 0.018$	$0.065\pm 0.018$	$0.067\pm 0.018$	$0.066\pm 0.018$
$n_s$	$0.9655\pm 0.0061$	$0.9448\pm 0.0070$	$0.9284\pm 0.0088$	${0.9191}_{-0.0078}^{+0.0076}$	${0.9161}_{-0.0072}^{+0.0081}$
$\ln \left({10}^{10}{A}_s\right)$	$3.089\pm 0.036$	$3.063\pm 0.035$	$3.023\pm 0.038$	$3.027\pm 0.037$	$3.028\pm 0.036$
$G_X/G_F$	–	${10}^8$
$m_s$	–			$[0.93,1.30]$	$1.27\pm 0.028$
$H_0$	$67.31\pm 0.95$	${62.2}_{-1.7}^{+2.0}$	${62.6}_{-1.8}^{+1.8}$	$59.56\pm 0.88$	${58.91}_{-0.79}^{+0.82}$

Table 3. Table 3 : Parameter constraints for the models under consideration, from the PlanckTT+lowP+BAO dataset. We either quote constraints in the form “mean ± plus-or-minus \pm 68% uncertainty”, or as 95% credible intervals (when not indicated, the lower limit should be understood to be zero). Units of m s subscript 𝑚 𝑠 m_{s} and H 0 subscript 𝐻 0 H_{0} are eV and km s -1 Mpc -1 , respectively.

Parameter	S$\Lambda$CDM
${\Omega }_b{h}^2$	$0.02197\pm 0.00021$
${\Omega }_c{h}^2$	${0.1144}_{-0.0015}^{+0.0016}$
$100{\theta }_{MC}$	${1.04332}_{-0.00063}^{+0.00090}$
$\tau$	$0.074\pm 0.018$
$n_s$	$0.9392\pm 0.0063$
$\ln \left({10}^{10}{A}_s\right)$	$3.038\pm 0.036$
$G_X/G_F$
$m_s$
$H_0$	$65.26\pm 0.68$

Table 4. Table 4 : Best-fit χ 2 superscript 𝜒 2 \chi^{2} values for the models under consideration, for the PlanckTT+lowP dataset.

Parameter	$\Lambda \mathrm{CDM}$	S$\Lambda \mathrm{CDM}$_GX0	S$\Lambda \mathrm{CDM}$	S$\Lambda \mathrm{CDM}$_Broad	S$\Lambda \mathrm{CDM}$_Narrow
${\chi }_{\min }^2$	$11265.1$	$11272.8$	$11269.0$	$11275.2$	$11277.6$

Equations

$$\mathcal{L}=g_X{\bar{\nu }}_s{\gamma }_{\mu }\frac{1}{2}(1-{\gamma }_5){\nu }_s{X}^{\mu },$$

$$\mathrm{i}\frac{d\rho }{dt}=[\mathsf{\Omega },\rho ]+C[\rho ],$$

$$\mathsf{\Omega }$$

$$G_X=\frac{\sqrt{2}}{8}\frac{g_X^2}{M_X^2}.$$

$${\Gamma }_t\simeq ⟨P({\nu }_{\alpha }\to {\nu }_s){⟩}_{\mathrm{coll}}{\Gamma }_X,$$

$${\Gamma }_X\simeq G_X^2{T}_{\nu }^5\frac{p}{⟨p⟩}\frac{n_s}{n_a}.$$

$$⟨P({\nu }_{\alpha }\to {\nu }_s){⟩}_{\mathrm{coll}}\simeq \frac{1}{2}{\sin }^2{\theta }_{\alpha s}.$$

$$(n_e,n_{\mu },n_{\tau },n_s{)}_{\mathrm{initial}}=(1,1,1,0),$$

$$(n_e,n_{\mu },n_{\tau },n_s{)}_{\mathrm{final}}=\left(\frac{3}{4},\frac{3}{4},\frac{3}{4},\frac{3}{4}\right),$$

$$G_X^2{T}_{\nu }^5\sim H(T_{\gamma }).$$

$$G_X\sim 10^{10}{G}_F,$$

$$\hat{L}[\delta f]=\hat{C}[\delta f],$$

$$\frac{\partial {\Psi }_i}{\partial \tau }+i\frac{q(\vec{k}\cdot \hat{n})}{ϵ}{\Psi }_i+\frac{d\ln f_0}{d\ln q}\left[\dot{\eta }-\frac{\dot{h}+6\dot{\eta }}{2}{\left(\hat{k}\cdot \hat{n}\right)}^2\right]=-{\Gamma }_{ij}{\Psi }_j,$$

$${\nu }_s\simeq \sin {\theta }_s{\nu }_1+\cos {\theta }_s{\nu }_4,$$

$${\Gamma }_{ij}=\left[\begin{array}{cccc}{\sin }^2{\theta }_s& 0& 0& \sin {\theta }_s\cos {\theta }_s\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ \sin {\theta }_s\cos {\theta }_s& 0& 0& {\cos }^2{\theta }_s\end{array}\right](3/2)(\zeta (3)/{\pi }^2)a{G}_X^2{T}_{\nu }^5.$$

$${\dot{\Psi }}_{i,0}=-\frac{qk}{ϵ}{\Psi }_{i,1}+\frac{1}{6}\dot{h}\frac{d\ln f_0}{d\ln q},$$

$${\dot{\Psi }}_{i,1}=\frac{qk}{3ϵ}\left({\Psi }_{i,0}-2{\Psi }_{i,2}\right),$$

$${\dot{\Psi }}_{i,2}=\frac{qk}{5ϵ}\left(2{\Psi }_{i,1}-3{\Psi }_{i,3}\right)-\left(\frac{1}{15}\dot{h}+\frac{2}{5}\dot{\eta }\right)\frac{d\ln f_0}{d\ln q}-{\Gamma }_{ij}{\Psi }_{j,2},$$

$$\displaystyle\dot{\Psi}_{i,\ell}=\frac{qk}{(2\ell+1)\epsilon}\Big{[}\ell\Psi_{i,(\ell-1)}-(\ell+1)\Psi_{i,(\ell+1)}\Big{]}-\Gamma_{ij}\Psi_{j,\ell}\quad(\ell\geq 3)\,,$$

Full text

1 Introduction

In recent years there has been a renewed interest towards light sterile neutrinos, suggested by different anomalies observed in short-baseline (SBL) neutrino experiments (see [1, 2, 3, 4] for recent reviews). In particular, these anomalies can be explained postulating a sterile neutrino with mass $m_{s}\simeq{\mathcal{O}}$ (1 eV) and active-sterile mixing angle $\theta_{s}\simeq 0.1$. For these values of the sterile neutrino parameters, the new states would be copiously produced in the Early Universe, resulting in a conflict with existing cosmological bounds on primordial radiation density and neutrino mass [5, 6, 7]. For this scenario to survive, a mechanism must be in place to suppress sterile neutrino abundance in the early universe: e.g., large primordial neutrino asymmetries [6, 8, 9], free primordial power-spectrum [10] or low reheating temperature [11]. Recently, a new mechanism has been proposed, which achieves such a suppression postulating secret interactions among sterile neutrinos, mediated by a massive gauge boson $X$, with $M_{X}\ll M_{W}$ [12, 13, 14] (see [15, 16] for the case of sterile neutrinos interacting with a light pseudoscalar). These secret interactions are described by the following Lagrangian

[TABLE]

where $g_{X}$ is the gauge coupling. Secret interactions would generate a large matter term in the sterile neutrino sector, that reduces the effective mixing angle, suppressing the active-sterile oscillations. Since the secret interactions are confined to the sterile sector, at the beginning this scenario seemed unconstrained. However, it was later realized that as the matter potential generated by the secret coupling declines as the Universe expands, sterile neutrinos would eventually encounter a resonance, when the matter potential becomes of the order of the neutrino vacuum oscillation frequency. This would allow for a sterile neutrino production through the combination of resonant Mikheyev-Smirnov-Wolfenstein effect [17, 18] and of non-resonant production via the secret collisions [19]. In this regard, it has been shown by some of us that for a coupling constant $g_{X}\gtrsim 10^{-2}$ and masses of the mediator $M_{X}\gtrsim{\mathcal{O}}$(10 MeV) the sterile neutrino production would occur before neutrino decoupling ($T\gtrsim$ 0.1–1 MeV) impacting the yield of light elements during Big Bang Nucleosynthesis (BBN) [20]. For smaller values of the mediator mass, BBN would be unaffected. However, in this case, sterile neutrinos would still be produced at $T\ll$ 0.1 MeV, when the matter potential becomes smaller than the vacuum oscillation term. Even assuming a negligible resonant production, sterile neutrinos would be copiously produced by the collisional term in the secret sector seeded by vacuum mixing, analogously to the Dodelson-Widrow mechanism acting for dark matter sterile neutrinos [21]. In [22] some of us have shown that this decoherent sterile neutrino production would quickly lead to equilibrium among active and sterile species, leading to a sizeable abundance of the latter in conflict with the cosmological neutrino mass bound. In addition, this mechanism reduces the effective number of neutrinos to $N_{\rm eff}\simeq 2.7$ at matter-radiation equality. For $M_{X}\gtrsim 0.1$ MeV, sterile neutrinos would be free-streaming before becoming non-relativistic and they would affect the structure formation at scales smaller than the free-streaming length. Conversely, for masses $M_{X}\lesssim 0.1$ MeV, as noticed in [22, 23], sterile neutrinos would be at the border between free-streaming and collisional regime at the photon decoupling, so one cannot naively apply the mass constraints as we did before. This range of the parameter space for the secret interactions is particularly interesting since it was previously shown [13, 14] to have potentially important consequences for the small scale structure of dark matter if the mediator $X$ couples also to dark matter. Furthermore, possible signatures of secret interactions in the observations of very-high-energy neutrinos by Icecube has been analyzed in [24, 25] (see also [26, 27, 28] for recent studies on the impact of secret interactions among active neutrinos on Icecube observations).

In the present work, we pursue a dedicated investigation of $M_{X}\lesssim 0.1$ MeV region obtaining constraints by the latest Planck data on the cosmic microwave background. The plan of this paper is as follows. In Sec. 2 we discuss the production mechanism of sterile neutrinos associated with secret interactions in the post-decoupling epoch and we present the existing cosmological bounds on this scenario. In Sec. 3 we present the results of our analysis and we draw our conclusions in Sec. 4.

2 Secret interaction framework

2.1 Sterile neutrino production

The 3+1 active-sterile neutrino mixing scenario involves 3 active families and a sterile species. Describing the neutrino system in terms of $4\times 4$ density matrices $\rho=\rho(p)$, the active-sterile flavour evolution is ruled by the kinetic equations [29]

[TABLE]

see [8] for a detailed treatment. The first term on the right-hand side of Eq. (2.1) describes the flavour oscillations Hamiltonian, given by

[TABLE]

where ${\sf M}^{2}$ = ${\mathcal{U}}^{\dagger}{\mathcal{M}}^{2}{\mathcal{U}}$ is the neutrino mass matrix in flavour basis, with ${\mathcal{U}}$ the active-sterile vacuum mixing matrix. The terms proportional to the Fermi constant $G_{F}$ in Eq. (2.2) are the standard matter effects in active neutrino oscillations, while the term proportional to $G_{X}$ represents the new matter secret potential. In particular, ${\sf E_{\ell}}$ is related to the energy density of $e^{-}$ and $e^{+}$ , ${\sf E_{\nu}}$ is the $\nu$-$\nu$ interaction term proportional to a primordial neutrino asymmetry (that here we assume negligible), while ${\sf E_{s}}$ is the energy density associated with $\nu_{s}$. The last term in the right-hand side of Eq. (2.1) is the collisional integral given by the sum of the standard ($\propto G_{F}^{2}$) and the secret one ($\propto G_{X}^{2}$). Since the flavour evolution typically occurs at neutrino temperature $T_{\nu}\ll M_{X}$ we can reduce the secret interaction to a contact form, with an effective strength

[TABLE]

The strong collisional effects produce a damping of the resonant transitions and would bring the system towards the flavour equilibrium among the different neutrino species with a production rate given by [30, 22]

[TABLE]

where $\langle P(\nu_{\alpha}\to\nu_{s})\rangle_{\rm coll}$ is the average probability of conversions among an active $\nu_{\alpha}$ and a sterile neutrino $\nu_{s}$ in a scattering time scale $(\Gamma_{X})^{-1}$, where the scattering rate is given by

[TABLE]

In Eq. (2.5) $\langle p\rangle\simeq 3.15T_{\nu}$ is the average-momentum for a thermal Fermi-Dirac distribution, and $n_{s}$ and $n_{a}$ the sterile and active neutrino abundance, respectively.

2.2 Cosmological bounds: state-of-the-art

Since the search for sterile neutrinos in laboratory experiments is presently open, it is important to use as many observations as possible to corner sterile neutrinos and in particular their production through secret interactions. In this context, cosmological observations represent a valid complementary tool to probe this scenario, being sensitive to the number of neutrinos, to their mass and to their free streaming characteristic.

In this Section we present the cosmological bounds obtained so far. For a coupling constant $g_{X}\gtrsim 10^{-2}$ and masses of the mediator $M_{X}\gtrsim{\mathcal{O}}$(10 MeV) the sterile neutrino production would occur before neutrino decoupling ($T\gtrsim$ 0.1–1 MeV). At this regard, in [20] it has been computed the sterile neutrino production relevant for BBN. The standard BBN dynamics is altered both by a larger value of $N_{\rm eff}$ and by the spectral distortion of $\nu_{e}$ when oscillations occur close to the neutrino decoupling. Using the present determination of deuterium primordial abundances, it was found that the 2H/H density ratio excludes much of the parameter space at 3$\sigma$, in particular masses $M_{X}\geq 40\leavevmode\nobreak\ \textrm{MeV}$ are excluded.

For smaller values of the mediator mass a large matter potential is generated suppressing the sterile neutrino production before the neutrino decoupling. With this choice of parameter ranges, BBN is left unchanged and gives no bound on the model. However, at lower temperatures when active-sterile oscillations are no longer matter suppressed, sterile neutrinos are still in a collisional regime, due to their secret self-interactions. The interplay between vacuum oscillations and collisions leads to a scattering-induced decoherent production of sterile neutrinos with a fast rate given in Eq. (2.4). At this regard, in [22] were neglected the resonant matter effects in the sterile neutrino production, reducing the average probability in Eq. (2.4) to a pure vacuum one, i.e.

[TABLE]

Taking as representative mixing angle $\sin^{2}2\theta_{es}\simeq 0.12$ [34], one would expect a sterile neutrino abundance, $n_{s}\simeq 0.06\,\ n_{a}$. This seemingly negligible population is enough to generate a large scattering rate in the post-decoupling epoch for a sufficiently large $G_{X}$. In particular for $G_{X}\gtrsim 10^{8}G_{F}$ the scattering rate at $T_{\gamma}\sim 10^{-2}$ MeV would be much larger than the Hubble rate $H$. This would lead a fast flavour equilibration between the three active and the sterile species, leading from an initial abundance

[TABLE]

to a final one:

[TABLE]

for all the parameters associated with eV sterile neutrino anomalies.

Soon after $\nu_{s}$ are produced via oscillation, active and sterile neutrinos have a shared grey-body distribution, namely a Fermi-Dirac function weighted by a factor 3/4 for each species. However, in the presence of strong secret interactions, these grey-body distributions will fastly evolve towards a Fermi-Dirac equilibrium function. The constant number density (or entropy) constraint implies that the temperature of this final spectrum is reduced by a factor $(3/4)^{1/3}$ with respect to the initial active neutrino temperature $T_{\nu}=(4/11)^{1/3}T_{\gamma}$. As a consequence of this effect, the total energy density stored in active and sterile neutrinos is reduced and the value of the effective number of neutrino species decreases down to $N_{\rm eff}\sim 2.7$ for relativistic neutrinos. A further slight reduction would occur at the matter radiation equality, i.e. for $T_{\gamma}\sim 0.7$ eV since eV sterile neutrinos would not be fully relativistic.

Secret interactions also affect the evolution of perturbations in the sterile neutrino fluid. In fact, if sterile states scatters via secret interactions, the free streaming regime is delayed until the scattering rate becomes smaller than the Hubble parameter. It means that if $G_{X}$ is large enough so that this condition holds at the non relativistic transition, sterile neutrinos would never have a free streaming phase, but always diffuse [22, 23]. One can obtain the smaller value of $G_{X}$ for which this occurs comparing the scattering rate with the Hubble rate $H$ at a temperature $3.15T_{\nu}\sim\langle p\rangle\sim m_{\rm s}$

[TABLE]

Writing the Hubble rate for $T_{\gamma}\sim$eV and using that $T_{\nu}=(4/11)^{1/3}(3/4)^{1/3}T_{\gamma}$ one obtains:

[TABLE]

which corresponds to $M_{X}\sim 10^{-1}$ MeV for $g_{X}\simeq 10^{-1}$.

In Ref. [22], sterile neutrino production through secret collisions has been examined also in view of the cosmological bounds on the sterile neutrino mass. Since for $G_{X}\lesssim 10^{10}G_{F}$ massive sterile neutrinos are free-streaming at the transition to the nonrelativistic regime, the peculiar effect of suppression of small-scale matter perturbations, induced by the presence of a light, free-streaming species, is preserved. In this case it is legitimate to expect that the mass bounds for non-interacting neutrinos coming from Cosmic Microwave Background (CMB) and Large Scale Structure (LSS) observations still apply also in the interacting case, as they basically rely on the effect of sterile neutrinos on the perturbation evolution. With this assumption, it was argued in Ref. [22] that an eV-mass sterile state, as suggested by the short baseline laboratory anomalies, was in tension, at least at the $2\sigma$ level with cosmological bounds available at the time, including those from the 2013 data release of the Planck satellite [35]. On the contrary, if the coupling is extremely large ($G_{X}>10^{10}G_{F}$), the free-streaming regime is reached only after the nonrelativistic transition, and the cosmological mass bounds possibly do not apply.

In this work, we aim to perform a dedicated, self-consistent analysis to derive observational bounds on secret contact interactions using the latest public Planck data on CMB anisotropies, as well as information from baryon acoustic oscillations (BAO). In order to do that, we do not rely on mass bounds obtained with the standard assumption of free-streaming neutrinos, but instead derive our own bounds, taking into account neutrino scattering via secret interactions. Moreover, we also take into account another effect induced by secret interactions, namely the increased density and pressure perturbations in the neutrino fluid, due to the fact that collisions cause power to flow towards the lower moments of the Boltzmann hierarchy. This produces a distinctive signature in the CMB anisotropy spectrum, as we shall see in the next section. Both these improvements are obtained by considering a collision term, proportional to the scattering rate, in the Boltzmann hierarchy for neutrinos, and by performing a fully consistent analysis in which both the sterile neutrino mass and the effective strength are treated as free parameters.

2.3 Secret interactions and cosmological perturbations

In this Section we describe the effect of interactions among neutrino species on the evolution of cosmological perturbations, that in turn determine the observational signatures on the CMB anisotropies and on large scale structures. Writing the sterile neutrino distribution function $f$ as the sum of a zero-th order part $f_{0}$ and a perturbation $\delta f\equiv f_{0}\Psi$, the latter evolves according to the Boltzmann equation:

[TABLE]

where $\hat{L}$ is the Liouville operator. The collision term $\hat{C}$ in the right-hand side takes into account the effect of secret interactions. In principle, the collision term is a complicated integral involving the matrix elements for the relevant processes; computing exactly the collision integral is a numerically demanding task, beyond the scope of our work (see e.g. Refs. [31, 32] for a detailed study of this topic). Fortunately it is enough, for the purpose of studying the evolution of cosmological perturbations, to resort to the so-called relaxation time approximation [33], in which the collision integral is taken to be $\hat{C}[\delta f]\simeq\delta f/\tau_{c}$, $\tau_{c}=\langle an\sigma v\rangle^{-1}$ being the mean conformal time between collisions. We can rewrite the Boltzmann equation in a more convenient way (we refer the reader to Ref. [36] for the notation):

[TABLE]

where the indices $i$ and $j$ label neutrino mass eigenstates, and summation over repeated indices should be understood. In the case under consideration, the scattering cross section $\sigma$ is of the order of $G_{X}^{2}T_{\nu}^{2}$, where $T_{\nu}=(3/11)^{1/3}T_{\gamma}$ is the common temperature of active and sterile neutrinos after flavour equilibration. Given that the neutrino number densities $n_{s}=n_{\nu}=(3/2)(\zeta(3)/\pi^{2})\,T_{\nu}^{3}$, we have that the collision rate $\Gamma=\tau_{c}^{-1}\sim a\,G_{X}^{2}T_{\nu}^{5}$. Comparing this with the conformal Hubble expansion rate $\mathcal{H}\equiv aH$, we can find the time at which collisions cease to be important and sterile neutrino start to behave as free-streaming particles.

Boltzmann codes like camb [37] evolve the perturbations in the distribution functions of the mass eigenstates. In order to obtain the scattering rates between mass eigenstates, those should be projected from the flavour basis through the mixing matrix. We shall assume that the sterile state is the superposition of the 1 and 4 mass eigenstates through the vacuum mixing angle $\theta_{s}$ as

[TABLE]

so that we are in the situation in which the mass eigenstates $\nu_{1}$ and $\nu_{4}$ interact with relative rates $\sin^{2}\theta_{s}$ and $\cos^{2}\theta_{s}$, while $\nu_{2}$ and $\nu_{3}$ are essentially free-streaming [23], and the scattering rate term becomes:

[TABLE]

It is possible to rewrite the Boltzmann equation for the mass eigenstates as an infinite hierarchy of multipoles [36]:

[TABLE]

where $\ell$ is the parameter of the Legendre expansion. We have set to zero the $\ell=0$ and $\ell=1$ terms of collision integral, in order to ensure particle number and momentum conservation. Thus the scattering directly affects the neutrino fluid from the shear onwards, and then propagates to the lower order moments. It is clear that, as long as the collision rate is much larger than the expansion rate, interacting neutrinos behave as perfect fluid.111In this case, the system of Eqs. (2.15) becomes stiff and a direct numerical integration is unfeasible. During this tight-coupling regime, we only evolve the $\ell=0,\,1$ moments of the distribution, using an approximate form for $\ell=2$ to close the hierarchy (see e.g. Ref. [38]).This means that shear and higher moments are exponentially suppressed, and the power in fluctuations is bound to the local monopole and dipole of the neutrino fluid. The net effect is that, at scales that are within the horizon during the interacting regime, density and pressure perturbations are enhanced with respect to the non-interacting case. This enhancement propagates to the photon fluid, and thus to CMB anisotropies, through the metric perturbations, as it can be clearly seen in Fig. 1, where we plot the temperature angular power spectrum (APS) (multiplied by an additional factor of $\ell^{2}$) for three models with an interacting sterile neutrino with $m_{s}=1\,\mathrm{eV}$ and different values of the coupling strength $G_{X}$. In all cases we take $N_{\mathrm{eff}}=2.7$, consistently with the expectation of flavour equilibration. The prediction for the case with $G_{X}\sim 10^{8}G_{F}\simeq 10^{3}\,\mathrm{GeV}^{-2}$ is practically identical to that of a $\Lambda\mathrm{CDM}$ extension with one non-interacting sterile neutrino and $N_{\mathrm{eff}}=2.7$. This means that there is a range of values around $G_{X}\sim 10^{8}G_{F}$ in which we still have a copious production of sterile neutrinos and flavour equilibration, but no direct effects of the interaction are visible on the APS (still, the effective number of relativistic degrees of freedom is $N_{\mathrm{eff}}=2.7$). Larger values of $G_{X}$ change the spectrum by increasing the power below a critical scale, related to the size of the horizon at the time at which neutrinos enter the free-streaming regime. For the parameter values used in the plot, we have that the comoving scale that enters the horizon at this time is $k\simeq 0.01\,\mathrm{Mpc}^{-1}$ ($k\simeq 0.03\,\mathrm{Mpc}^{-1}$), mapped to $\ell\simeq 130$ ($\ell\simeq 400$) for $G_{X}=10^{9}G_{F}$ ($G_{X}=10^{10}G_{F}$). In the following we will also consider BAO data to derive constraints on the parameters of the model. It is known that the inclusion of BAO measurements greatly improves constraints on neutrino masses, due to the breaking of geometrical degeneracies (most importantly, the one between the mass and $H_{0}$). As we shall see, BAO data also help in constraining the strength of secret interactions. Even if the effect of secret interactions is only seen at perturbation level, while BAO probe the background expansion, nevertheless they help to break parameter degeneracies that are present when only CMB are considered. In particular, the effect of a delayed neutrino decoupling can be canceled by increasing the total matter density, since the lesser amount of early-integrated Sachs-Wolfe effect will compensate the enhancement of perturbation power described above. We thus expect that constraining the matter density through BAO observations will result in tighter constraints on $G_{X}$.

3 Cosmological analysis

In this section we discuss our analysis of available cosmological data to constrain the coupling of the secret interaction. We first describe the method and data used, and then we present our results.

3.1 Method and data

We use the Boltzmann code camb [37], modified as described in the previous section, to follow the evolution of cosmological perturbations and compute the CMB anisotropy power spectra for given values of the cosmological parameters, including the secret coupling $G_{X}$ and the mass $m_{s}\equiv m_{4}$ of the (mostly) sterile neutrino. In order to derive Bayesian credible intervals for the parameters, we use the Markov Chain Monte Carlo (MCMC) engine CosmoMC [39] (interfaced with the modified camb). Our parameter space consists of the six parameters of the $\Lambda$CDM model, augmented by the parameters describing the sterile neutrino sector. The $\Lambda$CDM parameters are the baryon and cold dark matter densities $\omega_{b}\equiv\Omega_{b}h^{2}$ and $\omega_{c}\equiv\Omega_{c}h^{2}$, the angle $\theta$ subtended by the sound horizon at recombination, the optical depth to reionization $\tau$ and the logarithmic amplitude $\ln(10^{10}A_{s})$ and spectral index $n_{s}$ of the primordial spectrum of scalar fluctuations. The neutrino sector is instead described by the secret coupling $G_{X}$ and the sterile neutrino mass $m_{s}$. As explained in the previous section, we fix $N_{\mathrm{eff}}=2.7$, consistently with the assumption that all neutrino states (both active and sterile) have a common temperature $T_{\nu}=(3/11)^{1/3}T_{\gamma}$. This amounts to the assumption that $G_{X}\gtrsim 10^{8}G_{F}$. We also fix the active-sterile mixing angle to $\theta_{s}=0.1$ and the sum of the masses of the (mostly) active neutrinos to $0.06\,\mathrm{eV}$, equally shared among three mass eigenstates. We further assume flat spatial geometry and adiabatic initial conditions.

In our analysis, we always take flat, wide (in the sense that they are much larger than the expected posterior widths) priors for the six $\Lambda\mathrm{CDM}$ parameters. We also consider priors on $G_{X}$ and $m_{s}$ in order to model limiting cases of the scenario under consideration, to include additional pieces of experimental information, or simply to explore different regions of the parameter space. We start by performing a set of exploratory MCMC runs in which we assume a flat prior distributions in $\log_{10}\left[G_{X}\right]$ and $m_{s}$. The advantage of a logarithmic prior in $G_{X}$ is that it allows to explore several orders of magnitude in the parameter with equal probability per decade and thus to assess when the effect of secret interactions on the CMB APS becomes “large”, at least in comparison with the experimental sensitivity. However, a logarithmic prior gives more weight to small values of the parameter with respect to a flat prior, resulting in tighter bounds on the parameter itself. Moreover, it is an improper prior, since it does not integrate to a finite value if $G_{X}\geq 0$, and in order to give meaningful credible intervals an arbitrary, non-zero, lower bound on $G_{X}$ has to be assumed. For these reasons, we only use the results from this analysis to estimate the sensitivity of the data to $G_{X}$ and to gauge the initial step of the subsequent MCMC runs, that always use a flat prior on $G_{X}$.

The full model, in which the $\Lambda\mathrm{CDM}$ parameters as well as $G_{X}$ and $m_{s}$ are varied, is dubbed S$\Lambda\mathrm{CDM}$ (standing for “$\Lambda\mathrm{CDM}$ with secret interactions”). In this case, and unless otherwise stated, we take flat and wide priors also on $G_{X}$ and $m_{s}$. Note that, as explained above, we always have $N_{\mathrm{eff}}=2.7$. A limiting scenario is obtained by fixing $G_{X}$ to a very small value in our modified camb while keeping $N_{\mathrm{eff}}=2.7$, in order to reproduce the case in which $G_{X}$ is large enough for the flavour equilibration to happen, while still being small enough not to affect the evolution of cosmological perturbations. As we have mentioned, this approximately corresponds to $G_{X}\sim 10^{8}G_{F}$. Since, as noted in the previous section, this case is practically indistinguishable, as long as the evolution of cosmological perturbations in concerned, from a $\Lambda$CDM scenario with $N_{\mathrm{eff}}=2.7$ and $G_{X}=0$, we shall refer to this model as “S$\Lambda\mathrm{CDM}$_GX0”. Finally, we also consider prior on $m_{s}$ to model information from short baseline experiments. We refer to Ref. [3] in which the allowed $3\sigma$ (i.e., 99.73% CL) range for the squared mass difference $\Delta m^{2}_{41}=m_{4}^{2}-m_{1}^{2}$ that explains the SBL anomalies is $0.87\,\ \mathrm{eV}^{2}\leq\Delta m^{2}_{41}\leq 2.04\,\ \mathrm{eV}^{2}$. Not knowing the full shape of the probability density distribution for $m_{s}$, we decided to model it considering two “extreme” cases: in the first (“narrow $m_{s}$ prior”) we impose a gaussian prior $m_{s}=1.27\pm 0.03\,\ \mathrm{eV}$ (the width of the prior is chosen to match the 1$\sigma$ confidence interval for $\Delta m^{2}_{41}$ [3], assuming $m_{4}\gg m_{1}$), while in the second (“broad $m_{s}$ prior”) we impose a flat prior $0.93\,\ \mathrm{eV}\leq m_{s}\leq 1.43\,\ \mathrm{eV}$, corresponding to the $3\sigma$ interval reported above. Finally, we will often compare our results to those obtained in the framework of the standard $\Lambda\mathrm{CDM}$ model; for these, we refer to the values reported in the Planck 2015 parameters paper [40], and in this case it should be understood that $N_{\mathrm{eff}}=3.046$ [41, 42]. A list of the abbreviations used for the models considered in this paper, including a short description, can be found in Tab. 1.

Our data consists of the baseline Planck 2015 dataset (dubbed “PlanckTT+lowP” in the Planck papers), that includes temperature data in the range $2\leq\ell\leq 2500$, as well as the large-scale ($2\leq\ell\leq 30$) polarization (based on the measurements of the 70 GHz channel) [43]. The likelihood function associated to the data is computed using the code publicly released by the Planck collaboration222We acknowledge the use of the products available at the Planck Legacy Archive (http://www.cosmos.esa.int/web/planck/pla).. We marginalize over a number of nuisance parameters related to astrophysical foregrounds and instrumental uncertainties, as described in Ref. [43]. We also consider geometrical information coming from baryon acoustic oscillations; in particular we make use of the BAO results from the 6dF Galaxy Survey [44], from the BOSS DR11 LOWZ and CMASS samples [45], and from the Main Galaxy Sample of the Sloan Digital Sky Survey [46]. The extended dataset combining the Planck 2015 data with the BAO information will be denoted “PlanckTT+lowP+BAO”.

3.2 Results

We are now ready to present our results, summarized in Tabs. 2 and 3, where we show the Bayesian credible intervals for the parameters, for the various models and dataset combinations under consideration. As seen above, the presence of an interacting sterile neutrino impacts the cosmological observables in three ways:

• •

smaller $N_{\mathrm{eff}}$ due to flavour equilibration;

• •

larger density of (possibly) free-streaming species;

• •

reduced shear in the neutrino component of the cosmological fluid.

We start by considering the limit of small $G_{X}$ ($\sim 10^{8}G_{F}$), in which the third effect listed above is negligible, in order to disentangle the first two effects. Comparing the columns for $\Lambda\mathrm{CDM}$  and S$\Lambda\mathrm{CDM}$${\mathrm{GX0}}$, we note that there are considerable shifts in some parameters, in particular $H_{0}$ and $n_{s}$. The direction of the shifts is consistent with what we would expect given the well-known degeneracies of these parameters with both $N_{\mathrm{eff}}$ and the total density in light species. Looking at the $\chi^{2}$ values for the best-fit models, reported in Tab. 4, we find that S$\Lambda\mathrm{CDM}$${\mathrm{GX0}}$ performs worse in terms of goodness-of-fit, with a $\Delta\chi^{2}=7.7$ with respect to $\Lambda\mathrm{CDM}$. This is due to the low value of $N_{\mathrm{eff}}$ imposed by the flavour equilibration, while Planck data prefer a value closer to the standard expectation $N_{\mathrm{eff}}=3.046$. We also note that the mass of the sterile is constrained to be $m_{s}<0.82\,\ \mathrm{eV}$ at 95% CL, thus being in strong tension with the values suggested by the SBL anomalies.

The impact of secret interactions can be assessed by varying the coupling strength as a free parameter of the model. To this purpose we compare S$\Lambda\mathrm{CDM}$${\mathrm{GX0}}$ to S$\Lambda\mathrm{CDM}$, shown in columns 2 and 3 of Tab. 2. There are several points worth noticing: first of all, the constraints on the mass of the fourth eigenstate do not change with respect to the case of small $G_{X}$, thus remaining in tension with the preferred SBL solution. Secondly, secret interactions stronger than $G_{X}=2.8\times 10^{10}G_{F}$ are disfavoured, precluding the possibility of the collisional regime lasting after $z\sim\mathrm{few}\times 10^{3}$. Thus the scenario in the sterile neutrino starts to free stream only after recombination, is disfavoured. This is consistent with the fact that the bound on $m_{s}$ that we get is of the same order of magnitude as the ones obtained by the Planck collaboration in a minimal extension of the $\Lambda\mathrm{CDM}$ model. In that analysis, the effective mass $m_{s}^{\mathrm{eff}}\equiv 94.1\Omega_{\nu}h^{2}$ is used to parametrize the contribution of the sterile neutrino to the cosmological energy density. It is straightforward to see that, in our model, $m_{s}^{\mathrm{eff}}=(3/4)\,m_{s}$, so that in terms of the effective parameter the 95% upper bound for S$\Lambda\mathrm{CDM}$ reads $m_{s}^{\mathrm{eff}}<0.61\,\mathrm{eV}$. This should be compared with the result from the Planck collaboration for the same dataset (taken from the parameter tables available at the Planck Legacy Archive, see footnote 2), $m_{s}^{\mathrm{eff}}<0.88\,\mathrm{eV}$ at 95% CL. The two values cannot be directly compared, since the Planck analysis considers $N_{\mathrm{eff}}$ as a free parameter, with a prior $N_{\mathrm{eff}}\geq 3.046$, while in our analysis it is fixed to $N_{\mathrm{eff}}=2.7$. However, the tighter limit we find for $m_{s}^{\mathrm{eff}}$ is consistent with the lower value of $N_{\mathrm{eff}}$, given the direct correlation between the two parameters. In any case, the best-fit $\chi^{2}$ for S$\Lambda\mathrm{CDM}$ is still worse than $\Lambda\mathrm{CDM}$ ($\Delta\chi^{2}=3.9$) but yet better than S$\Lambda\mathrm{CDM}$${\mathrm{GX0}}$. When we also include information from BAO, we get tighter limits on $G_{X}$ and, especially, $m_{s}$, with 95% upper bounds of $1.97\times 10^{10}\,G_{F}$ and $0.29\,\mathrm{eV}$, respectively (see Tab. 3). In Fig. 2, we show the joint constraints and the marginalized one-dimensional posterior distributions for $G_{X}$ and $m_{s}$. For comparison, in the two-dimensional plot, we also indicate with a red star a model with $G_{X}=1.5\times 10^{10}G_{F}$ and $m_{s}=1\,\mathrm{eV}$, representative of the strong self-interacting scenario of Refs. [13, 23], that was argued to reconcile cosmological measurements and sterile neutrino interpretation of SBL anomalies (note that the other scenario considered in Ref. [23], with weak self-interactions, is not mapped by our analysis). In particular, a value $G_{X}\sim 10^{10}G_{F}$ roughly corresponds to the white band in the upper left part of Fig. 4 of Ref. [23] (at least down to the point where the 4-point approximation is valid, namely $M_{X}\sim 10^{-2}$ MeV and $g_{X}\sim 10^{-3}$) and the red star in that figure corresponds to a model with $G_{X}=1.5\times 10^{10}G_{F}$. We stress that, even if from this figure the scenario considered in Refs. [13, 23] seems to be excluded at the $\sim 3\sigma$ level for our most conservative choice of the dataset, i.e. PlanckTT+lowP, and even more strongly for PlanckTT+lowP+BAO, the actual statistical significance of the exclusion is somehow larger in both cases. A proper assessment should take into account that models with sterile secret interactions with $G_{X}>10^{8}\,G_{F}$ have $N_{\mathrm{eff}}=2.7$, a value that is itself mildly disfavoured with respect to the $\Lambda\mathrm{CDM}$ prediction of $N_{\mathrm{eff}}=3.046$. In the following paragraph, we will better quantify this statement, for the PlanckTT+lowP dataset, by comparing $\chi^{2}$ values between the best-fit models for $\Lambda\mathrm{CDM}$ and S$\Lambda\mathrm{CDM}$.

In order to better assess the (dis)agreement between Planck CMB observations and the SBL anomalies, also in the presence of secret interactions, we look at the fourth and fifth columns of Tab. 2, where we show the results for the cases in which we impose priors on $m_{s}$ that mimic the preferred SBL solution. For the S$\Lambda\mathrm{CDM}$_Broad model (column 4 of Tab. 2) we obtain almost the same constraint on the strength of the secret interaction we have obtained for the S$\Lambda\mathrm{CDM}$ scenario, in spite of the larger value of $m_{s}$ imposed by the prior. We notice however that the posterior distribution for $m_{s}$ is peaked in the lower edge of the prior, $m_{s}=0.93\,\mathrm{eV}$. In the S$\Lambda\mathrm{CDM}$_Narrow analysis, on the other hand, the larger a priori value of the sterile neutrino mass, $m_{s}\simeq 1.27$ eV, yields a looser constraint $G_{X}<4\times 10^{10}\,G_{F}$. For these two models, we see that the best-fit $\chi^{2}$ (computed on the PlanckTT+lowP dataset) is much worse than $\Lambda\mathrm{CDM}$: $\Delta\chi^{2}=11.1$ and $12.5$ for the “broad” and “narrow” priors, respectively. We argue that the inclusion of BAO information would make the tension even stronger, given the preference of that dataset for small values of the sterile neutrino mass. Finally, we notice how all models with non-standard interactions show a preference for values of $H_{0}$ even smaller than the one obtained in the framework of $\Lambda\mathrm{CDM}$ (see corresponding row of Tab. 2) further increasing the tension between CMB and direct estimates of the Hubble constant [47]; this is not captured by the $\chi^{2}$ figures reported above, that are computed on CMB data only. The increased tension is due in part to the low value of $N_{\mathrm{eff}}$, and, in the models with the SBL priors, by the large value of $m_{s}$; both effects, as per known degeneracies, push towards a smaller $H_{0}$.

4 Conclusions

In this work we have investigated, using Planck 2015 observations and a compilation of BAO data, the feasibility of cosmological models with sterile neutrinos, in addition to the three active states of the standard model of particle physics, with new, secret self-interactions mediated by a massive vector boson and confined in the sterile sector. This model has been proposed in order to alleviate the tension between the preferred solution of the SBL neutrino anomalies and cosmological observations, that disfavour a fourth fully thermalized neutrino species. Notably the effect of the new interactions would be to effectively dilute the density of both the active and sterile states (leading to an effective number of relativistic species $N_{\mathrm{eff}}=2.7$, more compatible with the Planck data). However, the mass of the sterile neutrino required to explain the SBL anomalies still appears to be too large with respect to the corresponding cosmological bounds. It was not clear a priori if and to what extent such bounds could be evaded thanks to the secret interactions that, if very strong, could significantly delay the onset of sterile neutrino free streaming.

Secret interactions also leave an imprint on the CMB spectra, by extending the collisional regime for the neutrino fluid. Using this effect, we have constrained the effective “Fermi constant” $G_{X}$ of the new interaction to be smaller than $2.8\times 10^{10}\,G_{F}$ at 95% CL from the Planck 2015 temperature and large-scale polarization data. This limit is improved to $2.0\times 10^{10}\,G_{F}$ at 95% CL when information from BAO are included. These results disfavour the range, corresponding to $G_{X}\gtrsim 10^{10}G_{F}$, in which the onset of sterile neutrino free streaming is delayed until after recombination, and cosmological mass bounds could be possibly evaded. In fact, our self-consistent analysis yields, at 95% CL, $m_{s}<0.82\,\mathrm{eV}$ and $m_{s}<0.29\,\mathrm{eV}$ from the Planck 2015 data alone and in combination with BAO, respectively, smaller than the value required to explain SBL anomalies, allowing to conclude that the tension between the SBL oscillation experiments and CMB observations still holds even in extended models with secret sterile neutrino interactions. Even disregarding BAO data, secret interactions with $G_{X}\gtrsim 10^{8}G_{F}$ are disfavoured with respect to standard $\Lambda\mathrm{CDM}$, by CMB data, due to their prediction of a low $N_{\mathrm{eff}}$. Moreover, CMB estimates of the Hubble constant $H_{0}$ in the secret interactions framework are smaller than their $\Lambda\mathrm{CDM}$ counterparts, thus increasing the tension with astrophysical measurements of the same quantity.

We summarize our findings in Fig. 3, where we show the parameter space excluded by our analysis in terms of the dimensionless coupling constant $g_{X}$ and the mediator mass $M_{X}$. The excluded region coincides with the whole region in which our assumptions hold and the approximations used are valid. The light and dark blue areas show the region excluded by our work. In particular, the light blue region corresponds to values of the interaction strength $G_{X}>2.9\times 10^{10}\,\mathrm{GeV}^{-2}$, thus larger than the $95\%$ upper limit on this parameter from Planck. In the dark blue region the range $10^{8}\,\mathrm{GeV}^{-2}<G_{X}<2.9\times 10^{10}\,\mathrm{GeV}^{-2}$ is still disfavoured by the neutrino mass bound. The red star is representative of the strong self-interacting scenario described in Refs. [13, 23]. The regions where the approximations used in our study become to break down are colored in gray. The horizontal band in light gray band on top indicates the non-perturbative regime ($g_{X}\gtrsim 0.1$) while the dark gray triangle on the bottom-left is where the interaction cannot be described as four-point interaction. This is obtained when the temperature at which $\nu_{s}$ are produced (approximated with Eq. (12) of [20]) is comparable or larger than the mediator mass $M_{X}$.

Our analysis has excluded the possibility of a single sterile neutrino with $\sim 1$ eV mass and $\sim 0.1$ mixing (as preferred by the SBL anomalies) with active neutrinos, having strong, four-fermion pointlike self-interactions. This is because it is not possible to hide the cosmological effects of such a large neutrino mass by means of a reduced free-streaming, without at the same time injecting too much extra power in the CMB angular power spectra. As it can be seen by comparing our Fig. 3 with Fig. 4 of Ref. [23] (please note that the quantities reported on the vertical axes of the two figures are related by $\alpha_{s}=g_{X}^{2}/4\pi$), the present analysis excludes the thin white band in the upper left part of Fig. 4 of Ref. [23] (dubbed there “strong self-interactions” region), that was regarded as being of particular interest as it could help explain the problems that arise at small scales in cold dark matter models of structure formation. On the other hand, the region in the lower part of Fig. 4 of Ref [23], corresponding to weak self-interactions, is not probed by our analysis. Even if a first exploratory study indicated this region as possible solution of the light sterile neutrino problem, recent refined calculations show that also this possibility seems to be ruled out (see Refs [25], and [48]), due to the X-mediated s-channel process leading to efficient sterile neutrino production. To conclude, we remark that our analysis assumes that the mass of the mediator is much larger than the temperatures relevant for the problem, and that self-interactions are perturbative. Moreover, we have not considered the case of two or more species of sterile neutrinos [49], although we argue that, in the case of complete thermalisation, they would be even more in tension with observations due to i) an even lower value of $N_{\mathrm{eff}}$, and ii) a larger density of interacting species, presumably resulting in a stronger bound on $G_{X}$.

Acknowledgments

We thank Basudeb Dasgupta and Joachim Kopp for reading our manuscript and for valuable comments on it. This paper is based on observations obtained with the satellite Planck (http://www.esa.int/Planck), an ESA science mission with instruments and contributions directly funded by ESA Member States, NASA, and Canada. We acknowledge support from ASI through ASI/INAF Agreement 2014- 024-R.1 for the Planck LFI Activity of Phase E2. We acknowledge the use of computing facilities at NERSC (USA). We acknowledge the use of the products available at the Planck Legacy Archive. The work of G.M. is supported by the Istituto Nazionale di Fisica Nucleare (INFN) through the “Theoretical Astroparticle Physics” project. The work of A.M. is supported by the Italian Ministero dell’Istruzione, Università e Ricerca (MIUR) and Istituto Nazionale di Fisica Nucleare (INFN) through the “Theoretical Astroparticle Physics” project. The work of N.S. has been supported by the German Research Foundation (DFG) under Grant Nos. KO 4820/1 1, FOR 2239, EXC-1098 (PRISMA) and by the European Research Council (ERC) under the European Union s Horizon 2020 research and innovation programme (grant agreement No. 637506, $\nu$Directions ).