1 Introduction

The unknown origin of neutrino masses and mixing together with the existence of the dark matter (DM) component of the Universe constitute our most significant experimental evidence for physics beyond the Standard Model (SM) and therefore the best windows to explore new physics. Neutrinos and DM also share an elusive nature with very weak interactions with the other SM particles. Indeed, neutrinos only participate in the weak interactions of the SM while all direct and indirect searches for DM interactions with the SM, other than gravity, are so far negative or inconclusive. A tantalising avenue of investigation is the possibility of a stronger connection between these two sectors. In this case, the best way to probe DM would be through the neutrino sector.

Several works have investigated the phenomenology of a dominant interaction between the neutrino and DM sectors and the possibility to probe DM through neutrinos both via its cosmological implications [1,2,3,4,5,6,7,8,9,10,11,12,13,14] as well as through indirect searches [14,15,16,17,18]. In the presence of this interaction, DM would no longer be collisionless, but able to scatter with neutrinos in the Early Universe, affecting matter density fluctuations. Moreover, the power spectrum would show a suppression at small scales [9, 10, 14] or even an oscillatory pattern [3,4,5, 8]. Indirect detection searches for DM annihilating to neutrinos in the galactic centre have also been performed at neutrino detectors and used to constrain DM–neutrino interactions [14,15,16]. The propagation of neutrinos through DM halos could be modified as well, leading to dips in supernova neutrino spectra due to resonant interactions with DM [19, 20], or affect the spectrum or isotropy of the high energy cosmic neutrinos observed by IceCube [21,22,23].

However, it is not straightforward to envision a scenario in which the neutrino–DM interactions dominate the DM phenomenology. Naively, gauge invariance dictates that the interactions of the left-handed (LH) SM neutrinos with DM will be equal to those of their charged lepton counterparts in the SU(2) doublets. In this case, the best windows to DM would instead be the charged leptons rather than the more elusive neutrinos.

In this work, we will investigate some gauge-invariant SM extensions that lead to sizeable neutrino–DM interactions, exploring if neutrino probes could dominate our sensitivity to the dark sector. This is actually a rather natural possibility. In fact, if DM does not participate in any of the SM gauge interactions, the natural expectation is that the strongest connection to DM will be via singlets of the SM gauge group. Indeed, if non-singlet fields were involved instead, the dimensionality of the operators linking the two sectors would have to increase in order to comply with gauge invariance. This reasoning leads to the three well-known SM portals to the dark sector: the “gauge boson portal” [24], the “Higgs portal” [25], and the “neutrino portal” [17, 26, 27]. The neutrino portal includes the addition of right-handed (RH) neutrinos \(N_R\), which makes this option particularly appealing in connection to the evidence of neutrino masses and mixing from neutrino oscillations.

Since the neutrino portal relies on the mixing between \(N_R\) and the light SM neutrinos to connect the neutrino and DM sectors, this mixing needs to be sizeable. In the “canonical” seesaw mechanism [28,29,30,31,32], the smallness of neutrino masses is explained through a large Majorana mass for \(N_R\) and the mixings are then similarly suppressed by the large scale. This option, which is rather natural from the point of view of neutrino masses, worsens the Higgs hierarchy problem [33]. An interesting alternative is to explain the smallness of neutrino masses via a symmetry argument instead [34,35,36,37,38,39]. Indeed, in models with an approximate lepton number (L) symmetry such as the linear [34, 35] or inverse [40] seesaw mechanisms, neutrino masses are suppressed by the small L-breaking parameters while light neutrino mixing with \(N_R\) is unsuppressed. In the present study, we will assume relatively large mixing angles noting that they can be compatible with neutrino masses, but we will not specify a concrete neutrino mass generation mechanism, since these small lepton number violating parameters, and hence light neutrino masses, will have no significant impact on the DM-related phenomenology.

We will consider fermionic DM and, more specifically, Dirac DM, which has the richest phenomenology when interacting with SM neutrinos. Indeed, the dominant term in the annihilation cross section to neutrinos is not velocity suppressed, and DM annihilations therefore lead to interesting signatures in indirect searches. Alternative scenarios with a Majorana, scalar, or vector DM candidate will lead to a velocity-dependent annihilation cross section to neutrinos [14]. While such possibilities are viable, they are difficult to probe experimentally at neutrino detectors. This is due to the fact that the DM velocity in the halo today is \(v_{\mathrm{{halo}}} = 10^{-3}c\) [41], which significantly reduces the annihilation rate to neutrinos.

The article is organised as follows. In Sect. 2, we summarise relevant experimental searches for DM and constraints coming from cosmology. In Sect. 3, we consider the simplest gauge-invariant scenario, in which DM is coupled directly to the full SM lepton doublet. In this case, as expected, the charged lepton probes tend to dominate the constraints on the DM parameter space. Further, in Sect. 4, we introduce the neutrino portal involving one new Dirac sterile neutrino N, which will communicate with the dark sector. We present two realisations of the neutrino portal, for scalar [42,43,44,45] and vector [46] interactions between the DM and N in Sects. 5 and 6, respectively. For both of them, we investigate the parameter space, demonstrating that current and future neutrino experiments have the dominant role in constraining it. Finally, we conclude in Sect. 7.

2 Constraints on interactions of DM with SM particles

In the next sections, we will explore the parameter space of different possible gauge-invariant ways to realise interactions of neutrinos with DM. For each realisation, we will investigate whether it is possible for these DM–neutrino interactions to play a dominant role in the DM phenomenology. In particular, we will address whether or not the DM relic abundance can be achieved via the DM–neutrino interactions and/or if indirect DM searches via its annihilation into neutrinos (probed at neutrino detectors) can be the dominant test of the model parameter space. We will use the observables presented in this section to place constraints on the parameter space of each scenario.

2.1 Indirect detection searches for DM annihilation to neutrinos

DM annihilating in high density regions such as the Milky Way can generate a significant monochromatic flux of neutrinos with energy \(E_\nu = m_\chi \), where \(m_\chi \) is the DM mass. This flux is proportional to the integral of the DM density squared along the line of sight and can be searched for in neutrino detectors such as Super-Kamiokande (SK) [47] or Borexino [48].

Several analyses that use neutrino detectors to probe the DM parameter space have been performed in the literature [14,15,16, 18, 49,50,51,52]. For small DM masses in the range 2–17 MeV, we can exploit the upper bound on the monochromatic antineutrino flux set by Borexino [53] and convert it to a conservative upper bound of \(\langle \sigma v_r\rangle \lesssim 10^{-22}\)–\(10^{-20}\) cm\(^3/\)s on the thermally averaged annihilation cross section \(\sigma \) multiplied by the relative velocity \(v_r\) of DM particles, as discussed in Ref. [14]. Likewise, between 10 and 200 MeV, SK can place an upper bound of \(\langle \sigma v_r\rangle \lesssim 10^{-25}\)–\(10^{-23}\) cm\(^3/\)s (depending on the DM mass) [14]. For DM with a mass between 1 GeV and 10 TeV annihilating in the galactic centre, the SK collaboration has performed a dedicated analysis and set an upper bound of \(\langle \sigma v_r\rangle \sim 10^{-24}\)–\(10^{-22}\) cm\(^3/\)s [50]. We will also consider the general upper bound on \(\langle \sigma v_r\rangle \) derived in Ref. [15] by calculating the cosmic diffuse neutrino signal from DM annihilations in all halos in the Universe and comparing it to the measured atmospheric neutrino background by Fréjus [54], AMANDA [55], and SK. This bound applies to \(m_\chi \) in the range between 100 MeV and 100 TeV and excludes \(\langle \sigma v_r\rangle \gtrsim 10^{-23}\)–\(10^{-21}\) cm\(^3/\)s (depending on \(m_\chi \)). As argued in Ref. [15], this bound could be improved by one or even two orders of magnitude with dedicated analyses by existing neutrino experiments such as SK.

The next generation experiment Hyper-Kamiokande (HK) [56] will be sensitive to approximately one order of magnitude smaller cross sections in this mass range. Indeed, with a 187 kton fiducial mass and an exposure time of 10 years, HK could probe the parameter space almost down to the relic density cross section (\(\langle \sigma v_r\rangle = 3\times 10^{-26}\) cm\(^3/\)s [57]). Possible improvements such as additional mass from a second tank together with Gd doping for background reduction would allow to probe beyond this value [51]. Similarly, the ESS\(\nu \)SB project [58] envisions a 500 kton fiducial water detector, MEMPHYS [59], that would have slightly better sensitivity than HK from the additional fiducial mass. Similarly, future DM and neutrino detectors such as DARWIN [60] and DUNE [61] will be able to further constrain the DM annihilation cross section to neutrinos. DARWIN will set stronger bounds for DM masses between 100 MeV and 1 GeV [62], while DUNE will be able to exclude thermal DM masses between 25 and 100 MeV [52].

Competitive constraints from DM annihilations in the Sun to neutrinos, or other SM particles that decay to neutrinos, have also been derived by neutrino detectors such as SK [63] and IceCube [64]. These exploit the higher DM concentration expected in the solar interior since it could capture DM particles from the halo via scatterings. In all the realisations under study we explore the connection between the DM and neutrino sectors with very suppressed interactions with the rest of the SM, in particular with quarks. Thus, in these scenarios, the Sun does not accrete DM particles effectively and the constraints from these searches do not apply.

2.2 Indirect detection searches for DM annihilation to charged leptons

DM interactions with charged leptons will always be present either at tree level, if DM couples to the full doublet, or at loop level in the neutrino portal scenarios. Therefore, we will take into account indirect detection searches for DM annihilations to charged leptons from the Fermi satellite [65], as well as from their imprint in the cosmic microwave background (CMB) as observed by Planck [66, 67].

2.3 Direct detection searches

DM will not couple directly to the quarks in any of the scenarios that we will discuss. Nevertheless, such couplings will arise at loop level in a similar way to the DM–charged lepton interactions. As we will see, bounds from direct detection experiments, such as XENON1T [68], are so stringent that they will still constrain the parameter space for large DM masses. Recently, direct detection of sub-GeV DM via scattering off electrons has gained significant attention [69,70,71,72]. We have also considered this process and found it to be sub-leading with respect to other relevant constraints.

2.4 Constraints from cosmology

If DM remains in thermal equilibrium with neutrinos during Big Bang nucleosynthesis (BBN), it can spoil its predictions [73, 74]. Similarly, the effective number of neutrinos, as constrained by CMB measurements, would be affected if DM remained in equilibrium after neutrinos decoupled from the photon plasma [75,76,77]. Thus, to avoid these two effects, we will not consider DM masses \(m_\chi < 10\) MeV. Moreover, DM–neutrino interactions can also have an effect in the formation of large scale structures (LSS) since, as DM particles scatter off neutrinos, they diffuse out and erase small scale perturbations. This effect leads to a suppression of the amount of small scale structures today. By comparing LSS predictions to observations, one can set an upper bound on the strength of the elastic scattering between DM and neutrinos [8, 78]. Nevertheless, for the models we are presenting in this work, the mixing between the sterile and SM neutrino suppresses the neutrino-DM elastic scattering and, consequently, its effect on LSS constrains regions of the parameter space already ruled out by CMB and BBN constraints [14].

3 Coupling to the full lepton doublet

In this section, we will study the simplest scenario, in which the neutrino–DM interaction arises from a direct coupling to the full SM SU(2) lepton doublet. In order to avoid specifying the nature of the mediator, we will adopt an effective field theory approach, simply adding a \(d=6\), 4-fermion interaction.

3.1 Model

Since the 4-fermion operator needs to involve two LH SM lepton doublets \(L_\alpha = (\nu _{\alpha L}, \ell _{\alpha L})^T\), \(\alpha = e,\mu ,\tau \), its Lorentz structure is fixed to be \(\overline{L_{\alpha }}\gamma ^{\mu }L_{\alpha }\). For definiteness we will assume a vector structure for the DM part. An axial coupling would instead lead to a velocity-suppressed DM annihilation cross section to neutrinos for both DM relic abundance and indirect searches. The cross section for DM annihilation to charged leptons would however have an additional term only suppressed by the lepton mass, and thus, it would tend to dominate over the annihilation cross section to neutrinos. Therefore, we will not consider this option in what follows.

The Lagrangian describing the neutrino–DM interaction is thus given by

(3.1)

where \(\chi \) is a Dirac fermion DM particle, and flavour diagonal couplings \(c_\alpha /\Lambda ^2\) between DM and the lepton doublets have been assumed in order to avoid new sources of flavour violation. For the effective description to be consistent we will require that \(\Lambda ^2/c_\alpha \gg m_{\chi }^2\). The simplest UV completion which leads to the \(d=6\) operator in Eq. (3.1) is via the exchange of a new heavy vector boson that couples both to \(\chi \) and \(L_\alpha \).

The Lagrangian in Eq. (3.1) implies that, in this naive gauge-invariant scenario, the coupling between the SM neutrinos and DM will be accompanied by a DM–charged lepton coupling of the same strength. Therefore, the strongest constraints on this model will typically come from indirect searches for DM annihilations to charged leptons. The DM relic abundance will also be set by its annihilation into leptons, either neutrinos or charged leptons, with the annihilation cross section given by

$$\begin{aligned} \langle \sigma v_r \rangle \approx \frac{c_{\alpha }^2 m_\chi ^2}{2\pi \Lambda ^4} \left( 1-\frac{m_{\alpha }^2}{4m_\chi ^2}\right) \sqrt{1 - \frac{m_\alpha ^2}{m_\chi ^2}}, \end{aligned}$$
(3.2)

where \(m_{\alpha }\) is the lepton mass for the different \(\alpha \) flavour.

3.2 Results

In Fig. 1, we show regions in the parameter space of the DM mass \(m_\chi \) and the new physics scale \(\Lambda \) excluded by different experiments. The blue line corresponds to the correct DM relic density \(\Omega _\mathrm {DM} h^2 = 0.1193 \pm 0.0009\) [66] obtained through the thermal freeze-out mechanism. This line has been computed with micrOMEGAs [79]. In the upper hatched region, the DM–lepton interaction would be too weak, leading to overclosure of the Universe (\(\Omega _\mathrm {DM} h^2>0.12\)). In the region below the blue line, the relic density is smaller than the observed DM abundance. If there are additional production mechanisms contributing to the DM density, this region is also viable.

Fig. 1
figure 1

Constraints on the DM mass \(m_{\chi }\) and the new physics scale \(\Lambda \). The upper and bottom-left panels correspond to couplings to only one of the lepton doublets (electron, muon, or tau), while the bottom-right panel corresponds to all three couplings being of equal strength. Along the blue line we recover the correct DM relic abundance from thermal freeze-out. The coloured shaded regions are excluded by different experiments, while the hatched areas correspond to prospective sensitivities of future experiments. The lower bound \(m_\chi \gtrsim 10\) MeV is set by observations of the CMB and BBN. See text for further details

Full size image

The constraints from indirect DM searches outlined in Sect. 2 are shown as different shaded regions. The light green (Planck [66, 67]) and orange (Fermi satellite [65]) regions correspond to the bounds from DM annihilation to charged leptons described in Sect. 2.2. The remaining shaded regions correspond to the constraints from DM annihilation to neutrinos as searched for in neutrino detectors and summarised in Sect. 2.1. In the upper-left panel of Fig. 1, we show in different colours the bounds coming from different neutrino experiments. The SK analyses [14, 50] are shown in red while the Borexino bounds [53] are displayed in yellow. The pink colour corresponds to the bounds from [49] obtained by combining the atmospheric neutrino data.Footnote 1 The dark red hatched region corresponds to prospective sensitivity of experiments on DM-electron scattering [71], while the blue, black, and green hatched regions correspond to prospects from different neutrino experiments as described in Sect. 2.1. In the following panels and in the rest of the paper we show all present indirect detection constraints from neutrino experiments in pink colour.

As can be seen in Fig. 1, the strongest constraints come from DM annihilation to charged leptons as probed by Fermi-LAT [65] for \(\chi \overline{\chi } \rightarrow \tau ^+\tau ^-\), \(\mu ^+\mu ^-\) and from Planck [66, 67] for \(\chi \overline{\chi } \rightarrow \ell ^+ \ell ^-\), \(\ell = e,\mu ,\tau \). The latter are in agreement with the results of Ref. [80], where, in particular, the dimension 6 operator given in Eq. (3.1) has been analysed. Indirect searches at neutrino detectors will always play a sub-leading role as long as annihilation to charged leptons is possible. Indeed, present constraints from DM annihilation to charged leptons are strong enough to rule out the entire allowed region of the parameter space that could lead to the correct DM relic density as long as the coupling to electrons is sizeable. However, if DM dominantly couples to the heavier lepton generations, allowed windows open up for \(m_\chi < m_\mu ~(m_\tau )\) (see the upper-right and bottom-left panels of Fig. 1). In this case, the DM relic density would be set by its annihilation to neutrinos, and the most relevant present constraints come from the results of SK and Borexino. The prospects for HK and DUNE would be very promising in these scenarios, allowing to probe most of the parameter space up to and beyond where the relic density is entirely explained by freeze-out based on neutrino interactions.

Regarding the constraints that could be set by the DM effects in the spectrum or isotropy of high energy cosmic neutrinos as observed by IceCube [21], these would lie in the region of the parameter space already excluded by the number of relativistic degrees of freedom in the early Universe [75,76,77].

From Fig. 1 it is clear that, as long as light DM couples to the electron doublet, this option for a neutrino–DM coupling is mostly ruled out by DM–electron interactions. However, if the DM coupling to \(L_e\) is negligible and DM dominantly couples to \(L_\mu \) and/or \(L_\tau \), the viable part of parameter space with \(m_\chi <m_\mu ~(m_\tau )\) can be probed by the neutrino experiments.

4 Coupling via the neutrino portal

Given the results of the previous section, we will now explore whether the neutrino portal option is able to lead to a rich DM-neutrino phenomenology without being in conflict with indirect searches involving charged leptons. The first necessary ingredient is to have sizeable mixing between the SM neutrinos and the new sterile neutrinos that will mediate the DM interaction. Therefore, the sterile-light neutrino mixing should not scale with the light neutrino masses, unlike in the canonical seesaw mechanism. Therefore, we will instead attribute the smallness of neutrino masses to an approximate lepton number (or \(B-L\)) symmetry rather than to a hierarchy of scales between the Dirac and Majorana masses. The new singlets will thus form pseudo-Dirac pairs since lepton number violation will necessarily be very small to account for the lightness of SM neutrinos. This is the case for instance in the popular “inverse” [40] and “linear” [34, 35] seesaw mechanisms based on such a symmetry.

As a simplifying assumption we will here consider the addition of only one (pseudo-)Dirac sterile neutrino that will serve as portal between the SM neutrinos and DM. Neglecting this small lepton number violation, the couplings between the SM and the new Dirac singlet neutrino are given by

(4.1)

where N is the Dirac sterile neutrino and \(\tilde{H} = i \sigma _2 H^*\), with H being the Higgs doublet.

Electroweak symmetry breaking gives rise to the neutrino Dirac mass term

$$\begin{aligned} \left( \overline{\nu _{\alpha L}},\, \overline{N_L}\right) M_\nu N_R + \mathrm {h.c.}, \end{aligned}$$
(4.2)

where \(M_\nu = (\lambda _\alpha v ,\, m_N)^T\) is the neutrino mass matrix and \(v = \langle H^0\rangle = 174\) GeV is the Higgs vacuum expectation value (vev). Diagonalising \(M_\nu M_\nu ^\dagger \) with a \(4\times 4\) unitary matrix U,

$$\begin{aligned} U^\dagger M_\nu M_\nu ^\dagger \, U = {{\,\mathrm{diag}\,}}\left( m^2_1, m^2_2, m^2_3, m^2_4\right) , \end{aligned}$$
(4.3)

we find the mass of the heavy neutrino to be

$$\begin{aligned} m_4 = \sqrt{m_N^2 + \sum _{\alpha } |\lambda _\alpha |^2 v^2}. \end{aligned}$$
(4.4)

As expected, the lepton number symmetry forbids light neutrino masses. In order to account for neutrino masses, small breaking of this symmetry via terms such as \(\mu \, \overline{N_L} N^c_L\) (inverse seesaw), or \(\lambda _\alpha ' \overline{L_\alpha } \tilde{H} N^c_L \) (linear seesaw) can be added. Since these small parameters would have negligible impact in the phenomenology of neutrino–DM interactions, we will not consider them in what follows.

The neutrino mixing matrix U, which relates LH flavour neutrino fields and the neutrino fields with definite masses as

$$\begin{aligned} \begin{pmatrix} \nu _{\alpha L} \\ N_L \end{pmatrix} = U \begin{pmatrix} \nu _{iL} \\ \nu _{4L} \end{pmatrix}, \quad \alpha = e,\mu ,\tau , \quad i = 1,2,3, \end{aligned}$$
(4.5)

has the form

$$\begin{aligned} U = \begin{pmatrix} U_{\alpha i} &{} U_{\alpha 4} \\ U_{s i} &{} U_{s4} \end{pmatrix}. \end{aligned}$$
(4.6)

The upper-left \(3\times 3\) block \(U_{\alpha i}\) would correspond to the Pontecorvo–Maki–Nakagawa–Sakata (PMNS) matrix once the small lepton number-breaking terms that induce neutrino masses are taken into account. Note that this matrix, being a \(3 \times 3\) sub-block of a larger unitary matrix will, in general, not be unitary. The upper-right \(3\times 1\) block \(U_{\alpha 4}\) describes the mixing between the active flavour neutrinos and the LH component of the heavy neutrino with mass \(m_4\). The last row of the matrix U specifies the admixture of each \(\nu _{jL}\), \(j=1,2,3,4\), in the LH sterile neutrino \(N_L\). As we will see in what follows, the DM-related phenomenology is driven by the mixing of active-heavy mixing matrix elements \(U_{\alpha 4}\). We will use the unitarity deviations of the PMNS matrix to constrain these mixings [81]. The mixing elements of interest are given by

$$\begin{aligned} U_{\alpha 4}= & {} \frac{\theta _\alpha }{\sqrt{1 + \sum _\alpha |\theta _\alpha |^2}} , \qquad U_{s4} = \frac{1}{\sqrt{1 + \sum _\alpha |\theta _\alpha |^2}}, \nonumber \\&\quad \sum _{i=1}^3 |U_{si}|^2 = \sum _{\alpha =e}^\tau |U_{\alpha 4}|^2, \end{aligned}$$
(4.7)

with \(\theta _\alpha = \lambda _\alpha v/m_N\). Note that, even though the SM neutrino masses have been neglected, the mixing with the extra singlet neutrino that will act as portal can still be sizeable. For definiteness we will fix the mixing to the different flavours to their \(1\sigma \) limit from Ref. [81], namely:

$$\begin{aligned} |\theta _e| = 0.031, \quad |\theta _{\mu }| = 0.011, \quad |\theta _{\tau }| = 0.044. \end{aligned}$$
(4.8)

In the following sections, we will explore two possible ways in which these Dirac neutrinos could couple to the dark sector and become portals between it and the SM neutrinos.

5 Neutrino portal with a scalar mediator

In this first example, we will assume that DM is composed of a new fermion, singlet under the SM gauge group, and that a new scalar mediates the Dirac neutrino–DM interactions.

5.1 Model

The Lagrangian of the model we will consider is given by

(5.1)

where \(\chi \) is a Dirac fermion DM candidate and S is a complex scalar. The fields \(\chi \) and S form the dark sector of the model (they are SM singlets), while N serves as a mediator between the dark sector and SM. The Lagrangian in Eq. (5.1) respects a global \(U(1)_L\) lepton number symmetry under which \(L_\alpha \), N, and \(S^*\) have the same charge and which protects the SM neutrino masses. Moreover, the Lagrangian respects a global \(U(1)_D\) dark symmetry, under which \(\chi \) and S have equal charges. This preserved symmetry ensures the stability of \(\chi \), if \(m_\chi < m_S\), where \(m_S^2 = \mu _S^2 + \lambda _{SH}v^2\) is the mass squared of the scalar S. For \(m_\chi > m_S\), the roles of \(\chi \) and S would change, and S would be a DM candidate. While this possibility is perfectly viable, it is more difficult to probe at neutrino detectors, as the DM annihilation cross section to neutrinos is velocity-suppressed. In what follows we assume that \(m_\chi < m_S\) and focus on fermionic DM.

This model was previously considered in Refs. [43, 45]. However, we will go beyond these works by performing a comprehensive analysis of the sensitivity of neutrino experiments to the parameter space of this model.

We will limit ourselves to the case in which DM is lighter than the heavy neutrino,Footnote 2 i.e., \(m_\chi < m_4\). This is the so-called direct annihilation regime [82], since DM annihilates through the mediator directly to SM particles. As intended, the only channel for DM annihilation at tree level is the one into light neutrinos. This process occurs via a diagram involving a t-channel exchange of the scalar mediator S. In the opposite regime, which is usually referred to as secluded [82], DM annihilates to heavy neutrinos, which subsequently decay. The phenomenology of this regime has been studied in Refs. [83,84,85,86].

Neglecting velocity-suppressed terms, we find the following thermally averaged cross section for DM annihilation to neutrinos:

$$\begin{aligned} \langle \sigma v_r\rangle\approx & {} \frac{y_L^4 }{32\pi }\left( \sum _{i = 1}^3 |U_{s i}|^2\right) ^2 \frac{m_\chi ^2}{\left( m_\chi ^2+m_S^2\right) ^2}\nonumber \\\approx & {} \frac{y_L^4 }{32\pi }\left( \sum _{\alpha =e,\mu ,\tau } \left| \theta _\alpha \right| ^2\right) ^2 \frac{m_\chi ^2}{\left( m_\chi ^2+m_S^2\right) ^2}. \end{aligned}$$
(5.2)

The product \(y_L \sqrt{\sum _\alpha |\theta _\alpha |^2}\) controls \(\langle \sigma v_r\rangle \) and, in order to allow for sufficient annihilation to reproduce the observed relic density, it cannot be too small. The value of the coupling \(y_L\) is limited by the requirement of perturbativity. We will restrict ourselves to \(y_L < 4\pi \). Since the coupling \(y_R\) does not enter Eq. (5.2), and thus, does not affect the tree-level DM–neutrino interactions, in what follows we set it to zero for simplicity. Regarding the mixing parameters \(\theta _\alpha \), the bounds on them depend on the mass of the heavy neutrino. For definiteness we will assume that the heavy neutrino has a mass above the electroweak scale. At this scale the bounds on heavy neutrino mixing derived in the global analysis of flavour and electroweak precision data performed in Ref. [81] apply. If smaller masses were instead considered, more stringent constraints from collider and beam-dump searches and, eventually, production in meson and beta decays could potentially apply [87] (see discussion in Sect. 6.3). In any case, all the observables relevant to DM phenomenology have a sub-leading dependence on \(m_4\). We also consider the case where the coupling \(\lambda _{SH} = 0\), ensuring the neutrino portal regime. In Refs. [43, 45], the radiative generation of the \(|S|^2 H^\dagger H\) operator was considered and its effects on \(m_S\) as well as on the invisible width of the Higgs boson were found to be negligible.

Fig. 2
figure 2

Thermally averaged annihilation cross section multiplied by the relative velocity for \(\chi \overline{\chi } \rightarrow \nu \overline{\nu }\). We have fixed \(m_S = 3m_\chi \), \(\theta _e = 0.031\), \(\theta _\mu = \theta _\tau = 0\), and varied \(y_L\) between 0.1 and \(4\pi \)

Full size image

In Fig. 2, we show the region of the parameter space for which the correct thermal relic abundance is obtained. This region spans DM masses up to 100 GeV for \(|\theta _e| = 0.031\), \(\theta _\mu = \theta _\tau = 0\), and \(y_L\) between 0.1 and \(4\pi \) while keeping \(m_S = 3m_\chi \) as a benchmark.

Annihilation of DM into charged lepton-antilepton pairs \(\ell ^+\ell ^-\) (\(\ell = e,\mu ,\tau \)) proceeds via the one-loop diagramsFootnote 3 shown in Fig. 3 (in unitary gauge).

The dominant contribution comes from the first and second diagrams, while the contribution from the last diagram is suppressed by the small Yukawa couplings of the charged leptons. The first diagram leads to the following effective operator:

$$\begin{aligned} \mathcal {L} \supset -a_{SW}\frac{g^2}{m_W^2} \overline{\chi } \gamma ^\mu P_R \chi \, \overline{\ell _\alpha } \gamma _\mu P_L \ell _\beta , \end{aligned}$$
(5.3)

where g is the weak coupling constant. Neglecting external momenta, the effective coupling \(a_{SW}\) is given by

$$\begin{aligned} a_{SW} = |U_{s4}|^2 U_{\alpha 4} U^*_{\beta 4} \frac{y_L^2}{(4\pi )^2} G\left( \frac{m_S^2}{m_4^2}\right) , \end{aligned}$$
(5.4)

where the loop function G(x) reads

$$\begin{aligned} G(x) = \frac{x - 1 - \log { x}}{4 \left( 1-x\right) ^2}. \end{aligned}$$
(5.5)

The second diagram in Fig. 3 leads to the following effective interaction of DM with the Z boson:

$$\begin{aligned} \mathcal {L} \supset - a_Z \frac{g}{\cos \theta _W} \overline{\chi } \gamma ^\mu P_R \chi Z_\mu , \end{aligned}$$
(5.6)

where \(\theta _W\) is the Weinberg angle and \(a_Z\) is the effective coupling, which in the limit of zero external momenta is given by

$$\begin{aligned} a_Z = |U_{s4}|^2 \left( 1 - |U_{s4}|^2\right) \frac{y_L^2}{(4\pi )^2} G\left( \frac{m_S^2}{m_4^2}\right) . \end{aligned}$$
(5.7)
Fig. 3
figure 3

One-loop diagrams (in unitary gauge) contributing to annihilation of DM into charged lepton-antilepton pairs \(\ell _\alpha \overline{\ell _\beta }\), \(\alpha ,\beta = e,\mu ,\tau \). The indices i and j run from 1 to 4

Full size image

These contributions have been also computed using a combination of packages: FeynRules [89, 90] to produce a model file, FeynArts [91] for generating the diagrams and FormCalc [92] for computing their numerical contributions. For numerical evaluation of the Passarino-Veltman functions we have used LoopTools [92]. We have also considered the limit of zero external momenta, which effectively corresponds to the limit of small DM and charged lepton masses, and confronted the analytical results obtained in this approximation using the package ANT [93] with the LoopTools results. For DM masses between 1 MeV and 100 GeV that we are interested in, the approximation works very well. The availability of analytical expressions allows for an easier exploration of the parameter space.

In Fig. 4, we present the cross sections for annihilation of DM into \(e^+e^-\), \(\mu ^+\mu ^-\), and \(\tau ^+\tau ^-\) for benchmark values of the model parameters. We fix \(m_S = 3m_\chi \), \(m_4 = 400\) GeV, \(y_L = 1\), \(\theta _e = 0.031\), and \(\theta _{\mu ,\tau } = 0\). As can be seen from the left panel, the annihilation cross sections to charged leptons are several orders of magnitude smaller than the cross section for DM annihilation into neutrinos. The difference in the cross sections becomes smaller when the DM mass approaches \(m_Z/2\), and the cross sections for \(\chi \overline{\chi } \rightarrow \ell ^+ \ell ^-\) exhibit a resonant behaviour due to the second diagram in Fig. 3. In the right panel, we show the indirect detection constraints from Planck [66, 67] and Fermi-LAT [65]. Note that those constraints assume a \(100\%\) annihilation rate into a single SM channel. Even for \(y_L = 4\pi \) the resulting annihilation cross sections into charged leptons are well below the experimental constraints. Thus, the considered realisation of the neutrino portal does provide an example of a gauge-invariant model in which the neutrino–DM interactions dominate DM phenomenology.

Fig. 4
figure 4

Thermally averaged annihilation cross section multiplied by the relative velocity for DM annihilation into \(e^+e^-\), \(\mu ^+\mu ^-\), and \(\tau ^+\tau ^-\). We have fixed \(m_S = 3m_\chi \), \(m_4 = 400\) GeV, \(y_L = 1\), \(\theta _e = 0.031\), and \(\theta _{\mu ,\tau } = 0\). The left panel provides comparison with \(\langle \sigma v_r\rangle \) for DM annihilation into neutrinos assuming the same set of model parameters. The right panel displays the indirect detection constraints coming from Planck and Fermi-LAT. The lower bound \(m_\chi \gtrsim 10\) MeV is set by observations of the CMB and BBN. See text for further details

Full size image

At one-loop level DM also interacts with quarks via diagrams involving Z and h, which are analogous to those in Fig. 3. The corresponding effective DM-nucleon spin-independent scattering cross section reads [45]

$$\begin{aligned} \sigma _n = \frac{\mu _n^2}{\pi } \frac{\left( Z f_p + \left( A-Z\right) f_n\right) ^2}{A^2}, \end{aligned}$$
(5.8)

where \(\mu _n\) is the reduced mass of the nucleon, A is the total number of nucleons in a nuclei, Z is the number of protons,

$$\begin{aligned} f_p = \left( 4 \sin ^2\theta _W - 1\right) \frac{G_F a_Z}{\sqrt{2}}, \qquad f_n = \frac{G_F a_Z}{\sqrt{2}}, \end{aligned}$$
(5.9)

with \(a_Z\) given in Eq. (5.7), and \(G_F\) being the Fermi constant. The radiative coupling of DM to the Higgs, \(\overline{\chi } \chi h\), would also give a contribution to direct detection searches. This contribution is however suppressed by the small quark Yukawa couplings. Direct detection of a SM singlet fermion DM candidate at one loop has been recently studied in detail in [94]. Moreover, an interesting example, which also provides radiative generation of neutrino masses, has been presented in [95].

The most stringent constraint on DM-nucleon spin-independent cross section for \(m_\chi \gtrsim 10\) GeV comes from XENON1T [68]. As we will see in the next subsection, this constraint is strong enough to probe the loop-suppressed scattering process if the value of the coupling \(y_L\) is sufficiently large. We have also considered DM scattering off electrons and found that the corresponding cross section is much smaller than the projected sensitivities of silicon, germanium, and xenon experiments derived in Ref. [71]. Thus, DM-electron scattering cannot provide an additional probe of the considered neutrino portal model.

5.2 Results

In this subsection, we explore the parameter space to find regions that satisfy all direct and indirect detection constraints and in which the DM phenomenology could be dominated by its interactions with SM neutrinos. We show our results in the \(m_\chi \)–\(m_S\) plane to determine the masses of the DM and the dark scalar that are presently allowed and could lead to the correct relic abundance (see Fig. 5).

Fig. 5
figure 5

Constraints on the DM mass \(m_\chi \) and the dark scalar mass \(m_S\). We have fixed \(\theta _e = 0.031\), \(\theta _{\mu ,\tau } = 0\); \(\theta _\mu = 0.011\), \(\theta _{e,\tau } = 0\); and \(\theta _\tau = 0.044\), \(\theta _{e,\mu } = 0\) (from top to bottom), considering \(y_L = 1\) and \(4\pi \). Along the blue line the DM relic density matches the observed value. The coloured shaded regions are excluded by different experiments, while the hatched areas correspond to prospective sensitivities of future experiments. The lower bound \(m_\chi \gtrsim 10\) MeV is set by observations of the CMB and BBN. See text for further details

Full size image

In Fig. 5, the triangular region \(m_S < m_\chi \) is forbidden by DM stability. Along the blue line(s) computed with micrOMEGAs,Footnote 4 the DM relic density matches the observed value \(\Omega _\mathrm {DM} h^2 = 0.1193 \pm 0.0009\) [66]. Above this line (the upper hatched region), the DM relic density is bigger than the measured value, i.e., DM overcloses the Universe. Below this line, the relic abundance would be smaller than the observed value. However, if there is an additional production mechanism, the relic abundance could also be compatible with this region.

As can be seen in the figure, indirect searches for annihilation to neutrinos, together with direct detection bounds by XENON1T for large DM masses, are the only probes that are presently constraining the allowed parameter space. The prospects to explore the remaining allowed regions through annihilation to neutrinos are very promising. In particular DUNE would be able to detect the neutrino signal in the range 25–100 MeV if the DM abundance is entirely due to this process.

In Fig. 6, we fix \(m_S\) to several representative values, namely \(m_S = 0.04\), 0.2, 1, and 5 GeV, and show the lines corresponding to the correct relic abundance in the \(m_\chi \)–\(y_L\) plane. These results have been obtained with micrOMEGAs. Small values of \(y_L\) are ruled out since they do not lead to efficient DM annihilation. As can be seen, a lighter dark scalar allows for smaller values of \(y_L\). For \(m_S \gtrsim 500\) MeV, the values of \(y_L \gtrsim 1\) are required to yield the observed relic density.

Overall, the cosmologically allowed parameter space of the model is already constrained by the current neutrino detectors as well as XENON1T.Footnote 5 Moreover, the next generation of neutrino experiments, in particular DUNE, will be able to probe thermal MeV fermion DM in the considered scenario.

6 Neutrino portal with a vector mediator

In this second example, we will again assume that DM is composed of a new Dirac fermion, this time coupled to a new massive vector boson. The Dirac singlet neutrino will also interact with this boson so as to provide the neutrino-DM interaction.

6.1 Model

The Lagrangian of the model is given by

(6.1)

where \(\chi \) is a Dirac fermion DM candidate, \(Z'\) is a new vector boson mediating the interaction between neutrinos and DM, and N is the Dirac sterile neutrino connecting the dark and visible sectors through its mixing with the active neutrinos. This Lagrangian could for instance describe a new \(U(1)'\) gauge symmetry spontaneously broken by the vev of a scalar SM singlet charged under it, that would induce masses for the \(Z'\) as well as for the heavy neutrino N and the DM. The particular mechanism is not relevant for the rest of the discussion and will not be elaborated further. We will also assume there is an additional conserved charge (e.g., a \(\mathbb {Z}_2\) symmetry) not shared between the neutrino and the DM that prevents their mixing. Note that in order to keep the Lagrangian in Eq. (6.1) anomaly free without introducing new fields, the simplest option is to couple the LH part of the Dirac sterile neutrino and the RH part of the DM to the new gauge boson with the same coupling \(g'\).

As in the previous scenario, we will assume that the DM mass \(m_{\chi }<m_4\) so that the dominant DM annihilation channel is to the three light SM neutrinos. This is a tree-level process and its cross section is given by

$$\begin{aligned} \left\langle \sigma v_r\right\rangle\approx & {} \frac{g'^4}{8\pi }\left( \sum _{i = 1}^3 |U_{s i}|^2\right) ^2 \frac{m_{\chi }^2}{(4m_{\chi }^2-m_{Z'}^2)^2} \nonumber \\\approx & {} \frac{g'^4}{8\pi }\left( \sum _{\alpha =e,\mu ,\tau } \left| \theta _\alpha \right| ^2\right) ^2 \frac{m_{\chi }^2}{(4m_{\chi }^2-m_{Z'}^2)^2}. \end{aligned}$$
(6.2)
Fig. 6
figure 6

Values of the DM mass \(m_\chi \) and the coupling \(y_L\) required to reproduce the observed relic abundance. We have fixed \(m_S = 0.04\), 0.2, 1, and 5 GeV, and have considered the representative case of \(\theta _e = 0.031\), while keeping \(\theta _{\mu ,\tau } = 0\). Along (above) the blue lines the DM relic density matches (is less than) the observed value. The lower bound \(m_\chi \gtrsim 10\) MeV is set by observations of the CMB and BBN

Full size image
Fig. 7
figure 7

Thermally averaged annihilation cross section multiplied by the relative velocity for \(\chi \overline{\chi } \rightarrow \nu \overline{\nu }\). We have fixed \(m_{Z'} = 3m_\chi \), \(\theta _e = 0.031\), \(\theta _\mu =\theta _\tau =0\), and varied \(g'\) between 0.1 and \(4\pi \)

Full size image

Note however that, for \(m_{Z'} \lesssim m_{\chi }\), the tree-level DM annihilation to a pair of \(Z'\) bosons is allowed. When this channel is open, it will dominate over the direct annihilation into neutrinos, since the latter is suppressed by neutrino mixing. This is the so-called secluded annihilation regime [82], which we do not consider in the present study.

Fig. 8
figure 8

One-loop diagrams contributing to the coupling of the \(Z'\) boson to charged leptons (left) and to kinetic and mass mixing between the \(Z'\) and Z bosons (right)

Full size image
Fig. 9
figure 9

Thermally averaged annihilation cross section multiplied by the relative velocity for DM annihilation into \(e^+e^-\), \(\mu ^+\mu ^-\), and \(\tau ^+\tau ^-\). We have fixed \(m_\chi {:}m_{Z_2}{:}m_4 = 1{:}3{:}6\), \(g' = 1\), \(\theta _e = 0.031\), and \(\theta _{\mu ,\tau } = 0\). The left panel provides comparison with \(\langle \sigma v_r\rangle \) for DM annihilation into neutrinos assuming the same set of model parameters. The right panel displays the indirect detection constraints coming from Planck and Fermi-LAT. The lower bound \(m_\chi \gtrsim 10\) MeV is set by observations of the CMB and BBN. See text for further details

Full size image

In this scenario, as can be seen from Fig. 7, the correct relic abundance can be obtained purely from annihilation to the SM neutrinos for values of the new gauge coupling \(g'\) between 0.1 and \(4\pi \), and DM masses in the 0.01–100 GeV range. In this figure, we have fixed \(m_{Z'} = 3 m_\chi \), \(|\theta _e| = 0.031\), and \(\theta _\mu = \theta _\tau =0 \) as benchmark values.

A direct coupling between the \(Z'\) boson and the charged leptons will also be induced through the loop diagrams in Fig. 8. Neglecting external momenta for the charged leptons, the effective vertex from the first loop diagram is given by

$$\begin{aligned} \mathcal {L}\supset - a_W g' \overline{\ell _\alpha } \gamma ^{\mu } P_L \ell _{\beta } Z'_{\mu }, \end{aligned}$$
(6.3)

where

$$\begin{aligned} a_W = |U_{s4}|^2 U_{\alpha 4}U_{\beta 4}^* \frac{g^2}{(4\pi )^2}\frac{m_4^2}{2m_W^2}. \end{aligned}$$
(6.4)

6.2 Mixing with the Z boson

Since the neutrino mass eigenstates have components that couple both to the Z and the \(Z'\), mixing between the two gauge bosons will be induced at loop level [24] through the second diagram in Fig. 8. The kinetic and mass mixings are described by the effective Lagrangian

$$\begin{aligned} \mathcal {L}_{Z'Z} = -\frac{\sin {\epsilon }}{2}Z'_{\mu \nu }Z^{\mu \nu }+\delta m^2Z'_{\mu }Z^{\mu }. \end{aligned}$$
(6.5)

Notice that these two terms could be present already at the Lagrangian level after gauge symmetry breaking. These would represent additional free parameters of the Lagrangian. However, these parameters do not contribute to the neutrino portal of interest here. Conversely, the neutrino mixing required for the neutrino portal does induce the Z–\(Z'\) mixing at the loop level. Barring fine-tuned cancellations between the allowed free parameters at the Lagrangian level and the loop-induced contributions from neutrino mixing, the minimum contribution present in our set-up will be the latter. We will therefore set the tree-level parameters to zero and require that the loop-induced contributions are below the present experimental constraints on Z–\(Z'\) mixing. We find the following results for the mixing parameters:

$$\begin{aligned} \delta m^2&= \frac{2}{(4\pi )^2}g'\frac{g}{\cos {\theta _W}}|U_{s4}|^2\left( 1-|U_{s4}|^2\right) m_4^2\, f_{1}, \end{aligned}$$
(6.6)
$$\begin{aligned} \sin {\epsilon }&= \frac{2}{(4\pi )^2}g'\frac{g}{\cos {\theta _W}}|U_{s4}|^2\left( 1-|U_{s4}|^2\right) f_{2}, \end{aligned}$$
(6.7)

where \(f_1\) and \(f_2\) are functions of \(x\equiv m_4^2/p^2\), namely,

$$\begin{aligned}&f_{1}(x) = \frac{1}{12}\bigg \lbrace 4x^2\left( 1-x^{-1}\right) ^3\nonumber \\&\quad \coth ^{-1}{(1-2x)}+2x-x^{-1}\log \left( {x}\right) -2\sqrt{x\left( 4-x^{-1}\right) ^3}\nonumber \\&\quad \arctan {\left( \left( 4x-1\right) ^{-1/2}\right) }\bigg \rbrace , \end{aligned}$$
(6.8)
$$\begin{aligned}&f_{2}(x) = -\frac{x^2}{6}\bigg \lbrace 4\left( 2x-3+x^{-2}\right) \coth ^{-1}{(1-2x)}\nonumber \\&\quad +4+x^{-2}\log {(x)} -2\sqrt{x^{-1}(4-x^{-1})}\left( 2+x^{-1}\right) \nonumber \\&\quad \arctan {\left( \left( 4x-1\right) ^{-1/2}\right) }\bigg \rbrace . \end{aligned}$$
(6.9)

For the purposes of this work \(p^2\sim m_{\chi }^2\), and thus, \(f_1\) and \(f_2\) will only depend on the ratio of the masses of the heavy neutrino and the DM particle. Following Ref. [96], we first diagonalise the kinetic term through a non-unitary transformation and then perform a rotation to diagonalise the mass term. The mass eigenstates \(Z_1\) and \(Z_2\) have masses given by

$$\begin{aligned} m_{Z_{1,2}}^2 = \frac{\sec ^2{\epsilon }}{2}\left( m_Z^2+m_{Z'}^2 -2 \delta m^2\sin {\epsilon } \mp \Delta \right) , \end{aligned}$$
(6.10)

where

$$\begin{aligned} \Delta&= {{\,\mathrm{sgn}\,}}\left( m_{Z'}^2-m_Z^2\left( 1-2\sin ^2{\epsilon }\right) -2 \delta m^2\sin {\epsilon }\right) \nonumber \\&\phantom {{}={}} \times \sqrt{m_Z^4+m_{Z'}^4+4 \delta m^4 -4 \left( m_Z^2+m_{Z'}^2\right) \delta m^2 \sin {\epsilon } -2 m_Z^2 m_{Z'}^2\left( 1-2\sin ^2{\epsilon }\right) }. \end{aligned}$$
(6.11)

From Eq. (6.10), one can easily verify that in the limit of small mass and kinetic mixing, i.e., \(\delta m^2 \rightarrow 0\) and \(\sin {\epsilon }\rightarrow 0\), the masses \(m_{Z_1}\rightarrow m_Z\) and \(m_{Z_2} \rightarrow m_{Z'}\). After the full diagonalisation, we can write the Z and \(Z'\) in terms of the mass eigenstates \(Z_1\) and \(Z_2\) as follows:

$$\begin{aligned} Z_{\mu }&=\left( \cos {\xi }-\tan {\epsilon } \sin {\xi }\right) Z_{1\mu } -\left( \sin {\xi }+\tan {\epsilon }\cos {\xi }\right) Z_{2\mu }, \end{aligned}$$
(6.12)
$$\begin{aligned} Z'_{\mu }&=\sec {\epsilon }\left( \sin {\xi }\, Z_{1\mu }+\cos {\xi }\, Z_{2\mu }\right) , \end{aligned}$$
(6.13)

where \(\xi \) is the angle related to the mass diagonalisation, which is defined through

$$\begin{aligned} \tan {\left( 2\xi \right) } = \frac{2\cos {\epsilon } \left( m_Z^2 \sin {\epsilon }-\delta m^2\right) }{m_{Z'}^2-m_Z^2\left( 1-2\sin ^2{\epsilon }\right) -2\delta m^2\sin {\epsilon }}. \end{aligned}$$
(6.14)

The two angles \(\xi \) and \(\epsilon \) will control the phenomenology associated to the Z-\(Z'\) mixing and consequently, the possible \(Z'\) couplings to fermions.

The loop-induced kinetic mixing parameter \(\epsilon \) depends solely on the ratio \(x \approx m_4^2/m_\chi ^2\), providing the coupling \(g'\) and the element \(U_{s4}\) of the neutrino mixing matrix are fixed (see Eqs. (6.7) and (6.9)), and increases with it. Fixing \(|\theta _e| = 0.031\) and \(\theta _{\mu ,\tau } = 0\), we find that for \(x = 4\), which is the lowest value preventing the \(\chi \overline{\chi } \rightarrow \nu _i\overline{\nu _4}\), \(i=1,2,3\), channels, and \(g' = 1~(4\pi )\), the mixing parameter \(|\sin \epsilon |\) is of order of \(10^{-6}~(10^{-5})\). For values of x as large as \(10^4\) and \(g' = 1~(4\pi )\), the value of \(|\sin \epsilon |\) does not exceed approximately \(10^{-5}~(10^{-4})\).

Generally, these values can be probed in beam dump and fixed target experiments searching for visible decay products (electrons and muons) of the \(Z_2\) boson with mass between approximately 1 MeV and 1 GeV (see, e.g., [97, 98]). However, in the considered model the \(Z_2\) decays mostly invisibly, either to a pair of the SM neutrinos or, if it is heavy enough, to a pair of DM particles, while its decays to charged leptons are suppressed. Thus, the bounds from fixed target experiments will not apply in this case. The supernova constraints cover nearly the same \(Z_2\) masses, but a different range of \(\epsilon \sim 10^{-10}\)–\(10^{-7}\) [97], which thus are also avoided. For larger \(Z_2\) masses, up to 100 GeV, collider experiments place the best constraints on \(\epsilon \sim 10^{-4}\)–\(10^{-3}\) (see, e.g., Ref. [98]). There exist also collider searchers for \(Z_2\) decaying invisibly, which constrain \(\epsilon \lesssim 10^{-3}\) for \(m_{Z_2} < 8\) GeV [99]. These collider constraints are above the values of the loop-induced kinetic mixing parameter in our model. Finally, the much weaker constraint from the invisible \(Z_1\) width, \(\epsilon \lesssim 0.03\) [100], is also evaded.

Together with the first diagram in Fig. 8, the size of \(\xi \) and \(\epsilon \) will determine how relevant the DM annihilation to a pair of charged leptons is. We find that the tree-level annihilation to neutrinos dominates over that to charged leptons. In Fig. 9, we show a particular example of this behaviour for \(m_4 = 2m_{Z_2}\), \(m_{Z_2} = 3m_\chi \), \(g'=1\), \(|\theta _e| = 0.031\), and \(\theta _\mu = \theta _\tau = 0\). It is clear from this figure that the annihilation to charged leptons is unconstrained by current experimental searches. Note that the Planck and Fermi-LAT constraints shown in the right panel of Fig. 9 assume a \(100\%\) annihilation rate into a single SM channel.

6.3 Results

The allowed regions of the parameter space in the \(m_\chi \)–\(m_{Z_2}\) plane that satisfy cosmological, indirect and direct detection constraints for this model are presented in Fig. 10 for \(g'=1\) and \(4\pi \), setting \(\theta _\alpha \ne 0\) one at a time and keeping two other mixing angles fixed to zero. For definiteness, in the figure we set \(m_4 = 2 m_{Z_2}\). Notice that this choice is not relevant for the interaction between the SM neutrinos and DM and only plays a role in the loop-induced processes that are sub-dominant. Nevertheless, if the \(Z_2\) originates from a new \(U(1)'\) gauge group, its mass \(m_{Z_2}\), as well as that of the Dirac neutrino \(m_4\), are generated after the breaking of the symmetry. Thus, the natural expectation is that \(m_4\) is not much heavier than \(m_{Z_2}\) as long as the new gauge coupling \(g'\) is \(\mathcal {O}(1)\). Hence, unlike for the scalar example, it is not appropriate to set \(m_4\) to a value above the electroweak scale while exploring (sub-)GeV \(Z_2\) boson masses.

Fig. 10
figure 10

Constraints on the DM mass \(m_\chi \) and \(m_{Z_2}\). Along the blue lines, computed with micrOMEGAs, the DM relic density matches the observed value. The coloured shaded regions are excluded by different experiments. The lower bound \(m_\chi \gtrsim 10\) MeV is set by observations of the CMB and BBN. See text for further details

Full size image

Below the electroweak scale constraints on the neutrino mixing parameters \(\theta _\alpha \) are a priori much more stringent [87]. However, in the model under investigation the heavy neutrino decays mostly invisibly to either a SM neutrino and the \(Z_2\) (if \(m_4 > m_{Z_2}\)), or a SM neutrino and a pair of the DM particles (if \(m_4 < m_{Z_2}\)), assuming \(g' \gtrsim 1\). This implies that the existing collider and beam dump constraintsFootnote 6 should be rescaled with the corresponding branching ratios and become even weaker than the non-unitarity constraints imposed previously for the scalar realisation. The bounds from peak searches in leptonic decays of pions and kaons will however apply, since they rely entirely on the kinematics of a two-body decay. Thus, the non-unitarity constraints actually dominate down to \(m_4 \approx m_K \approx 0.5\) GeV, where \(m_K\) is the kaon mass. In the region \(m_4 \sim 0.01\)–0.4 GeV, the bounds on \(U_{e4}\) and \(U_{\mu 4}\) from peak searches are very stringent. We do not display them explicitly in Fig. 10, because they are \(m_4\)-dependent, while all the constraints shown in the figures have an extremely sub-leading dependence on \(m_4\), as outlined above. Thus, Fig. 10 is to be interpreted as generally valid for any neutrino mass \(m_4 > m_K\).

The blue line was calculated with micrOMEGAs and represents the DM and vector boson masses that will produce the correct relic abundance in a thermal scenario, while the masses in the upper hatched area would generate too much DM. A key difference with respect to the previous model is that here the DM annihilation cross section to neutrinos proceeds via an s-channel and thus is enhanced for \(m_{Z_2} \sim 2 m_\chi \), as can be seen from Eq. (6.2). This explains the second branch of the blue line below the resonant condition in the panels with \(g'=1\). A line where the relic abundance can be obtained below \(m_{Z_2} = 2 m_\chi \) also occurs for \(g'=4\pi \) but, since the cross section is larger, the relic abundance is achieved for \(m_{\chi } > 100\) GeV, which is ruled out by XENON1T. This resonant effect also explains the shape of the indirect detection constraints which follow the same trend.

Similar to the previous model in Sect. 5, the direct detection constraints from XENON1T become relevant at large DM masses for \(g'=4\pi \). However, even for values of the gauge coupling this large, we have checked that direct detection constraints from the elastic DM scattering off electrons are negligible.

The complementarity between cosmological observables, DM, and neutrino experiments allows us to set very strong bounds on the DM and \(Z_{2}\) masses for this particular realisation, ruling out significant portions of the parameter space. There are still allowed regions for larger values of the gauge coupling consistent with a thermal DM candidate that yields the observed DM relic abundance. However, future neutrino experiments such as DUNE will be able to probe down to the value for which the correct relic abundance is obtained in some parts of the parameter space.

Fig. 11
figure 11

Values of the DM mass \(m_\chi \) and the coupling \(g'\) required to reproduce the observed relic abundance. We have fixed \(m_{Z_2} = 0.04\), 0.2, 1, and 5 GeV, and have considered the representative case of \(\theta _e = 0.031\), while keeping \(\theta _{\mu ,\tau } = 0\). Along (above) the blue lines the DM relic density matches (is less than) the observed value. We do not consider \(m_{\chi } > m_{Z_2}\) to ensure the neutrino portal regime. The lower bound \(m_\chi \gtrsim 10\) MeV is set by observations of the CMB and BBN

Full size image

It is worth noticing that the sensitivity of present and future neutrino detectors to DM annihilations into neutrinos is largely independent of the flavour to which the sterile neutrino dominantly couples. Indeed, regardless of the original flavour composition produced by the DM annihilations, neutrino oscillations will tend to populate all flavours with similar fractions when the flux arrives to the detector. The main differences between the three rows in Fig. 10 are due to the different magnitude of the mixing allowed to the different flavours, with more stringent constraints applying for the mixing with muon neutrinos.

Finally, in Fig. 11, we fix \(m_{Z_2}\) to several values, namely, \(m_{Z_2} = 0.04\), 0.2, 1, and 5 GeV, and show the lines corresponding to the correct relic abundance in the \(m_\chi \)–\(g'\) plane. These results were obtained using micrOMEGAs. Small values of \(g'\) are ruled out except for DM masses in the proximity of the resonance, i.e., when \(m_{\chi } \approx m_{Z_2}/2\). As can be seen from this figure, a lighter dark vector boson allows for smaller values of \(g'\). For \(m_{Z_2} \gtrsim 1\) GeV, values of \(g' \gtrsim 1\) are required to yield the observed relic density, except for the resonance region. The dip towards \(m_\chi \approx m_{Z_2}\) corresponds to opening of new DM annihilation channels at tree level.

7 Conclusions

Despite the tremendous improvement over the last years in the sensitivity of direct, indirect and collider searches for dark matter, its discovery still eludes us. An interesting possibility is that its interactions with SM particles happen dominantly with the neutrino sector. This option would not only explain our failure to detect any DM interactions (except gravitational) so far, it would also connect our two present experimental signals of physics beyond the SM. Indeed, a rich phenomenology that would stem from the connection of these two sectors has been explored and discussed in the literature. SU(2) gauge invariance would naively dictate that neutrinos share all their interactions with their charged lepton counterparts, which are much easier to detect. We have therefore explored whether a dominant neutrino–DM interaction is allowed in simple gauge-invariant models without conflicting with searches through charged leptons.

We first explored the simplest scenario, in which DM couples to the full lepton doublet. We verified that, as long as the DM is heavier than the charged lepton(s) it couples to, the bounds from DM annihilation to charged leptons preclude DM–neutrino couplings sizeable enough to be probed, even ruling out all of the parameter space that would not lead to overclosure of the Universe. Alternatively, if DM couples to \(\tau \) (\(\mu \)) and is lighter than the charged lepton, its phenomenology is dominated by the interaction with neutrinos. This region is constrained by present neutrino detectors and will be fully probed for certain DM masses by future experiments.

We have then explored the option of the neutrino portal to DM and showed, as an example, two specific realisations with scalar and vector couplings, respectively. In the neutrino portal DM couples directly to new heavy neutrinos. Indeed, their singlet nature makes them natural candidates to probe the dark sector since they are allowed to interact with it via relevant or marginal operators. These right-handed neutrinos are also a natural addition to the SM particle content so as to account for the evidence for neutrino masses and mixings. The mixing between the SM neutrinos and the new singlets will induce DM–neutrino interactions at tree level, but DM-charged lepton couplings only at loop level.

In the two realisations explored we find that it is indeed possible for neutrino detectors to place the most stringent and competitive bounds through searches for DM annihilations to neutrinos. Present searches at Super-Kamiokande, Fréjus, or Borexino are ruling out large areas of the parameter space. Interestingly, future projects such as Hyper-Kamiokande, MEMPHYS, DARWIN, or DUNE will be able to probe the cross section very close and beyond the value required to explain the DM abundance solely by annihilation to SM neutrinos. These new searches will effectively cover most of the parameter space, probing if the right-handed singlet fermions that can explain the origin of neutrino masses also represent our best window to the discovery of the dark matter sector.