1 Introduction

There are compelling reasons to believe that neutrino masses are the first manifestation of a new physics (NP) scale, which can be identified with the mass of the heavy mediator(s) that generates neutrino masses. Under this assumption, neutrino masses and mixings provide information on a combination of the mediator mass and couplings. Although an upper limit of the new physics scale can be derived by requiring the new couplings to be perturbative, generically, no lower bound results from this constraint. As it is well known, a very high NP scale leads to a hierarchy problem [1, 2], and to vacuum instability issues [3, 4], both of which can be avoided if the NP scale is not much higher than the electro-weak scale. The possibility that the neutrino mass mediators are light enough to be produced in accelerators and in atmospheric showers is therefore worth exploring.

The most popular extension of the Standard Model (SM) realizing neutrino masses is the Type-I seesaw model [5,6,7,8], with at least two heavy Majorana singlets \(N_j\):

$$\begin{aligned} {\mathcal {L}_N \supset \sum _j i \bar{N}_j /\!\!\!\partial N_j - \left( Y_{\alpha j} \bar{L}_\alpha {{\tilde{\Phi }}} N_j {+} {m_{N_j} \over 2} \bar{N}_j N_j^c \right) {+} \mathrm{h.c.} ,}\nonumber \\ \end{aligned}$$
(1)

where \({\tilde{\Phi }} \equiv i\sigma _2 \Phi ^*\) is the complex conjugate of the Higgs field \(\Phi \), \(L_\alpha \) is the SM lepton doublet with flavor \(\alpha \), Y is a generic Yukawa matrix, and \(m_{N_j}\) is the Majorana mass of the singlet \(N_j\). After spontaneous electro-weak symmetry breaking, the heavy Majorana states mix with the standard neutrinos resulting in a spectrum of three light states with masses \(m_\nu \propto m_N^{-1} (Y v)^2\), and two or more heavy neutral leptons (HNL) with masses \(\propto m_N\). In this model, all massive neutrino states are admixtures of the standard neutrinos and the singlet states, as dictated by the mixing matrix \(U_{\alpha j}\) (which diagonalizes the mass Lagrangian of the whole system). The phenomenology of the HNL depends crucially on their mass and their mixing to the charged leptons. In fact, it is through this mixing that the heavy singlets can be produced either through charged-current (CC) or neutral-current (NC) processes and also how they decay back to SM particles. For simplicity, in this work we will adopt the simplified notation \(|U_{\alpha }| \equiv U_{\alpha j}\), assuming that only one of these states is kinematically accessible for our purposes.

HNLs have been extensively searched for in laboratory experiments, using mainly two types of signatures: either displaced vertices from the decay of the HNL, or through precise determination of the decay product kinematics in meson decays (see, e.g., Refs. [9,10,11,12] for reviews of available constraints in the MeV–GeV mass range, or Refs. [13,14,15,16,17,18,19,20,21,22,23] for future prospects to improve over current bounds). In Ref. [24] we studied in detail the search for long-lived particles produced in the atmosphere, which would then decay inside large-volume neutrino detectors. We derived bounds on HNL for masses above the kaon mass, where laboratory limits are weaker. We found, however, that the limits from atmospheric searches are not competitive with laboratory constraints in this mass range.

In this letter we focus on the lighter mass region instead, where the HNL can be produced in kaon and pion decays for which the atmospheric flux is significantly larger. The most stringent bounds for HNL below the kaon mass come from peak searches [25, 26] and displaced HNL decay vertex searches [27,28,29]. However, for masses below the kaon and pion mass the HNL becomes very long-lived: although the value of its lifetime in the rest frame (\(\tau \)) depends strongly on its mass and on its mixing with the light states, in the minimal model described above it ranges between \( c\tau \sim (10^{-4} - 50 ) \times |U_\alpha |^{-2}\) (km), for \( m_N\) between 40 MeV and 400 MeV. This makes atmospheric neutrino detectors a well-suited setup to search for their decay products.

In this work, we use the framework developed in Ref. [24] to extract limits from the multi-GeV muon and electron neutrino data samples observed at Super-Kamiokande (SK). The search of atmospheric HNL in a similar mass range has been considered before for SK [30, 31] and IceCube [32]. Our analysis significantly improves the methodology of these earlier studies. In addition, extensions of the minimal model of Eq. (1) that decouple production and decay have been considered in recent works, particularly in relation with the LSND/MiniBoone anomaly [33,34,35,36,37,38]. In Refs. [39,40,41,42] an extension that includes a dipole interaction of the heavy singlets provides a new radiative decay channel, which dominates the HNL decay and significantly reduces its lifetime. Therefore, we will present our limits not only in the context of the minimal HNL model of Eq. (1) (that is, on the plane \(|U_\alpha |\) versus \(m_N\)), but also on the plane \(\mathrm{Br}(K/\pi \rightarrow N)\) versus \(c\tau \), which is phenomenologically motivated. As we will see, this can be useful in order to constrain non-minimal scenarios with uncorrelated production and decay, such as for example the dipole extension of Ref. [42].

2 HNL production mechanisms

The leading mechanism for HNL production is through the decays of mesons produced in the atmosphere. The computation is explained in detail in Ref. [24], and here we summarize it for convenience. Defining \(\ell \) as the distance traveled between the production point of the HNL to its point of entry in the detector, the HNL production profile (in a differential of distance \(d\ell \)) can be obtained as:

$$\begin{aligned} {d \Pi ^{\mathrm{decay}}_N \over d E d \cos \theta d \ell } = \sum _P \int _{E_P^{\mathrm{min}}}^{E_P^{\mathrm{max}}} {d E_P\over \gamma _P \beta _P c \tau _P} {d\Phi _P\over dE_P d\cos \theta } {d n \over d E }, \end{aligned}$$
(2)

where dn/dE stands for the differential distribution of the HNL energies while \(\Phi _P\) is the parent meson flux and, in this work, \(P= K^{\pm }, \pi ^\pm \). We have used the Matrix Cascade Equation (MCEq) Monte Carlo software [43, 44] to compute the fluxes for the parent mesons in the atmosphere, with theFootnote 1 SYBILL-2.3 hadronic interaction model [45], the Hillas-Gaisser cosmic-ray model [46] and the NRLMSISE-00 atmospheric model [47]. In Eq. (2), \(\beta _P\) and \(\gamma _P\) are the parent boost factors while \(\tau _P\) is its lifetime in the rest frame. For the dominant two-body decay \(K^\pm , \pi ^\pm \rightarrow N l^\pm _\alpha \) (denoted generally as \(P \rightarrow N Y\)), the differential distribution dn/dE reads

$$\begin{aligned} {d n ( E_P, E)\over d E }= & {} {1\over \Gamma _P} {d\Gamma (P\rightarrow N Y) \over dE} \nonumber \\= & {} {\mathrm{Br}(P\rightarrow N Y) \over \Gamma (P\rightarrow N Y)} {d\Gamma (P\rightarrow N Y) \over dE}, \end{aligned}$$
(3)

where Br stands for branching ratio, and

$$\begin{aligned} {1 \over \Gamma (P\rightarrow N Y)} {d\Gamma (P\rightarrow N Y) \over dE} = {1\over p_P \sqrt{\lambda (1,y_N^2, y_Y^2)}}, \end{aligned}$$
(4)

with \(y_i \equiv {m_i\over m_P}\) and

$$\begin{aligned} \lambda (a,b,c) \equiv a^2+ b^2 + c^2 -2 a b -2 a c -2 b c. \end{aligned}$$
(5)

Finally, the kinematical limits for \(E_P \gg m_P\) are:

$$\begin{aligned} E_P^{\mathrm{max(min)}} \equiv {2 E\over (1+ y_A^2-y_Y^2) \mp \sqrt{\lambda (1,y_A^2, y_Y^2)}} . \end{aligned}$$
(6)

Figure 1 shows a representative example of the HNL production from kaon and pion decays in the atmosphere, and compares it to the result obtained for heavier parent mesons and \(\tau \) leptons in our previous work [24]. As can be seen from this figure, the production profile grows by several orders of magnitude when the mass of the HNL allows for it to be produced from the decays of lighter resonances, due to their more abundant fluxes in the atmosphere.

Fig. 1
figure 1

HNL production rate (for \(m_N=0.1\) GeV) per energy, solid angle and path length, at a height of 15.4 km and \(\cos \theta = 0.9\), for different parent particles P. In this work, we consider only \(\pi \) and K decays. For comparison we also show the corresponding result for two-body decays of heavier mesons \((D, D_s)\) and leptonic (three-body) \(\tau \) decays from Ref. [24]. For all lines shown here, effects due to the mass of the charged lepton produced together with the HNL have been neglected for simplicity

Full size image

A second contribution to the flux comes from the HNL production in NC scattering of standard atmospheric neutrinos as they pass through the Earth’s matter (for instance in \( \nu \mathcal {N} \rightarrow N X\), where \(\mathcal {N}\) stands for a nucleon and X is a hadronic shower). This contribution can be estimated as follows. Assuming a flux of standard neutrinos \(\phi _\nu \), the HNL produced at a distance r from SK is:

$$\begin{aligned} {d \Pi ^{\mathrm{int}}_N \over d \cos \theta d E d r} {=} \int _{E_\nu ^{\mathrm{min}}}^{E_\nu ^{\mathrm{max}}} \!\!\!dE_\nu {d \phi _\nu (r)\over d \cos \theta d E_\nu } {d \sigma _{\nu \rightarrow N}(E_\nu ) \over d E} N_A\rho _{\oplus },\nonumber \\ \end{aligned}$$
(7)

where \(\sigma _{\nu \rightarrow N}\) is the HNL production cross section (which can be estimated as the standard \(\nu \) NC scattering cross section multiplied by the corresponding mixing \(|U_\alpha |^2\)), while \(N_A\) is the Avogadro number and \(\rho _{\oplus }\) is the Earth density at distance r. This assumes that both standard neutrinos and the resulting HNL are highly boosted and therefore the zenith angle of the HNL is roughly the same as that of the incoming neutrino. The upper and lower limits in the integral correspond to the kinematically allowed ranges, where a SM neutrino scattering on a nucleon at rest can produce a HNL with energy E.

3 HNL decays in Super-Kamiokande

The flux that arrives to the detector for HNL produced in meson decays, \(\Phi ^{\mathrm{decay}}_N\), is obtained integrating over all LLPs produced at different distances \(\ell \), weighted by their corresponding survival probabilities, as

$$\begin{aligned} {d \Phi ^{\mathrm{decay}}_N\over d\cos \theta dE} = \int _0^{\ell _{\mathrm{max}}} d \ell {d \Pi ^{\mathrm{decay}}_{N}\over d\cos \theta dE d\ell } ~e^{-{\ell \over L_{\mathrm{decay}}}}, \end{aligned}$$
(8)

where \(L_{\mathrm{decay}}\) is the decay length of the HNL in the laboratory frame. The maximum distance \( \ell _{\mathrm{max}} \equiv \ell (h_{\mathrm{max}},\theta )\) is a simple function of the maximum height of the atmosphere where cosmic showers are produced, \(h_{\mathrm{max}} \simeq 80\) km, and the zenith angle. Analogously, the flux entering the detector for HNL produced in the interaction of SM neutrinos in the Earth, \(\Phi ^{\mathrm{int}}_N\) is:

$$\begin{aligned} {d \Phi ^{\mathrm{int}}_N \over d\cos \theta d E} = \int _0^{R_{\mathrm{max}}(\theta )} dr {d \Pi _N \over d\cos \theta d E d r} e^{-r/L_{\mathrm{decay}}}, \end{aligned}$$
(9)

where \(R_{\mathrm{max}}(\theta ) = 2 R_{\oplus } \cos \theta \) is the maximum distance traveled through the Earth for trajectories with a zenith angle \(\theta \), \(R_\oplus \) being the Earth’s radius.

The total number of HNL decays inside the detector, within a given time window \(\Delta T\), energies in the interval \([E, E + dE]\) and zenith \([\cos \theta , \cos \theta + d\cos \theta ]\), can be computed as

$$\begin{aligned} \frac{d N}{dE d\cos \theta } = \Delta T A^{\mathrm{eff}}_{\mathrm{decay}}(E,\cos \theta ) {d \Phi _N\over dE d\cos \theta }, \end{aligned}$$
(10)

where \( \Phi _N \equiv \Phi ^{\mathrm{int }}_N + \Phi ^{\mathrm{dec}}_N\), and \( A^{\mathrm{eff}}\) is an effective area which accounts for the probability that a decay takes place inside the detector. This area can be estimated integrating the surface of the detector normal to the flux direction, weighted by the decay probability of the N inside the detector:

$$\begin{aligned} A^{\mathrm{eff}}_{\mathrm{decay}}(E,\cos \theta ) = \int dS_{\perp } \left[ 1- e^{- \Delta \ell _{det} (\cos \theta )\over L_{\mathrm{decay}}(E) } \right] \, . \end{aligned}$$
(11)

Here \(\Delta \ell _{det}\) is the length of the segment of the HNL trajectory that cuts into the detector (for explicit expressions we refer the reader to Ref. [24]). A cylindrical geometry for SK with height of 40 m and radius of 20 m is assumed.

The two contributions to the total number of events (coming from meson decays and from SM neutrino interactions in the Earth) have a very different angular dependence: while the flux from decays is expected to be larger from above, those from interactions come obviously from below. However, we have checked that the final contribution to the number of events coming from HNL produced in meson decays is several orders of magnitude larger than the one obtained from interactions of SM neutrinos on the Earth. The ratio between the two contributions decreases for larger \(c\tau \) but, within the range \(c\tau \in [1, 10^4]\) km, it lies within the range \([10^6, 10^2]\). We can therefore safely neglect the contribution from neutrino interactions in the Earth and, in the rest of this work, we will only consider HNL production from meson decays.

In order to derive our limits, we use the data samples as well as the expected neutrino background prediction from Fig. 5 of Ref. [48]. We use both the \(\mu \)- and e-like fully-contained multi-GeV events, adding the single- and multi-ring samples together (labeled as “multi-GeV” and “multi-ring” in Fig. 5 in Ref. [48], respectively). In the case of e-like events, we add both the \(\nu _e\) and \({\bar{\nu }}_e\) samples as well. In Ref. [48] data are binned in zenith angle, while information on the energy of the events is not publicly available. Therefore, in this work the data is binned only in \(\cos \theta \). As for the neutrino energies considered, we integrate over all energies between 1 GeV and 90 GeV, as this is the range corresponding to the fully-contained sample (see Fig. 6 in Ref. [48]). We think this is conservative, as SK may be sensitive to events outside this range.

The number of events observed in a given sample will also depend on detector efficiencies and reconstruction effects, which should be included in the form of migration matrices giving the relationship between true and reconstructed variables. Such information is however not publicly available for SK. Thus, here we make the simplifying assumption that the efficiencies are independent on the neutrino energy and, in particular, we take a flat detection efficiency \(\epsilon ^{\alpha } = 0.75\) both for \(\mu \) and e-like, in line with the values quoted in Ref. [48] for the multi-GeV \(\nu _e\) event sample. While a priori some loss of efficiency could be expected at high energies (mainly due to a reduction in the containment of the events), we think that the impact on our results would be small since we expect our sensitivities to come mostly from the events at low energies. This is because the heavy neutrino flux produced in the amtosphere will follow a very steep power law, peaked at low energies as shown in Fig. 1, where we expect the SK efficiencies to be best. We also assume that the angular reconstruction is much better than the width of the bins in zenith angle, so migration between different bins in \(\cos \theta \) can be neglected. If this assumption were to be relaxed, our sensitivities would probably be worsened. This is so because in the case of very long-lived particles (as is the case for HNL with masses below 500 MeV, and weakly coupled to the SM) their very long lifetimes lead to a higher sensitivity for specific angular bins, where the distance traveled by the HNL is comparable to (or smaller than) its lifetime in the laboratory frame. While here we have used only publicly available information, a detailed evaluation by the experimental collaboration is needed to validate our results. Finally, we should also point out that we expect an increase in sensitivity if the analysis were to be carried out using both energy and angular information. This is beyond the scope of this work.

The number of events in the i-th bin in \(\cos \theta \) in a given sample can therefore be computed as:

$$\begin{aligned}&N^{\mathrm{SK, \alpha -like}}_i =\mathrm{Br}(N \rightarrow \alpha \mathrm{-like}) \nonumber \\&\quad \times \int _{1~\mathrm{GeV}}^{90~\mathrm{GeV}} dE \epsilon ^{\alpha }(E) \int _{\cos \theta _{i}^{\mathrm{min}}}^{\cos \theta _i^{\mathrm{max}}} \! d\cos \theta \frac{dN}{dE d\cos \theta } \, , \nonumber \\ \end{aligned}$$
(12)

where \(\cos \theta _{i}^{\mathrm{min}}\) and \(\cos \theta _{i}^{\mathrm{max}}\) are the lower and upper limits of the bin. Here, \(\mathrm{Br}\)(\(\alpha \)-like) stands for the total branching ratio for all decay channels including muons, electrons or photons in the final state, depending on the sample (\(\alpha \)) considered. For example, in the case of \(\mu \)-like we consider only those decay channels including one or more muons in the final state. In the case of e-like events we require that no muons are present, but we also include decay channels with photons as these are easily mis-identified with electrons at SK (such as \(N \rightarrow \nu \pi ^0\), since the \(\pi ^0\) decays promptly to two photons).

4 Results

In this section we derive limits on HNL production by performing a \(\chi ^2\) fit to the SK data. A Poissonian \(\chi ^2\) function has been used:

$$\begin{aligned} \chi ^2 =2 \sum _{i, \alpha } N_i^{\alpha }+ B_i^{\alpha } - n_i^{\alpha }+n_i^{\alpha } \log \left( {n_i^{\alpha }\over N_i^{\alpha }+ B_i^{\alpha }}\right) , \end{aligned}$$
(13)

where the sum runs over the angular bins, and \(\alpha =\{e\)-like, \(\mu \)-like\(\}\). Here, \(n_i^\alpha \) stands for the data observed in each bin while \(N_i^\alpha \) is the predicted number of signal events and \(B_i^\alpha \) is the background prediction, which includes the predicted number of atmospheric neutrino events in the SM. In our \(\chi ^2\) calculations, we consider separately the e- and \(\mu \)-like samples and we add their two contributions to the total \(\chi ^2\). Therefore, our limits will be derived using 20 degrees of freedom, corresponding to the total number of bins in \(\cos \theta \). In our calculations we find, however, that the sensitivity is largely dominated by the e-like contribution since the size of the branching ratio \(\mathrm{Br}\)(e-like) is much larger than \(\mathrm Br\)(\(\mu \)-like) in the mass range considered.

First we show in Fig. 2 the results on the plane \(\mathrm{Br}(K/\pi \rightarrow N)\times \mathrm{Br}(N\rightarrow \mathrm{visible})\), where \(\mathrm{Br}(N \rightarrow \mathrm{visible})\) accounts for the probability that the HNL decays visibly in the detector. Our results are presented as a function of \(c\tau \), assuming no correlation between the production and decay mechanisms for the HNL. It is interesting to note that SK can outperform the powerful displaced-decay search limits from beam dumps such as those in Refs. [27,28,29], for models where the decay does not involve two charged tracks, since all laboratory searches request this condition to reduce background contamination. In the extended model of Ref. [42] this is precisely the case, since the HNL decay is dominated by the radiative \(N\rightarrow \nu \gamma \) decay via the dipole interaction, therefore the stringent limits from PS191 [27, 28] and the recent T2K limits [29] do not apply. The shaded purple region in Fig. 2 shows the range where the MiniBooNE anomaly could be explained, for a HNL with \(m_N=260\) MeV, extracted from Ref. [42]. In this case, the HNL would only be produced in kaon decays and, unfortunately, our limits from \(K^\pm \) decay (pink lines) fall short to probe this region. However, SK would be sensitive to non-minimal models with HNL produced in \(\pi ^\pm \) decays with a similar value of the production BR and lifetime and, obviously, larger neutrino experiments such as DeepCore or Hyper-Kamiokande could significantly improve over these constraints.

Fig. 2
figure 2

Limits at 90% confidence level (CL), on the plane \(\mathrm{Br}(K/\pi \rightarrow N) \times \mathrm{Br}(N \rightarrow \mathrm{visible})\) versus \(c\tau \) for \(m_N= 0.1\) GeV (solid) and \(m_N=0.25\) GeV (dashed), for HNL produced in \(\pi \) (dark blue) or K (pink) decays. Our limits have been obtained assuming uncorrelated production and decay mechanisms for the HNL. The shaded light blue region corresponds to the (correlated) values obtained for the minimal model, varying the mixings within their presently allowed constraints as explained in the text. The purple area corresponds to the \(1\sigma \) allowed region where a HNL with a dipole interaction could explain the MiniBooNE anomaly, according to Ref. [42], for \(m_N=260\) MeV

Full size image
Fig. 3
figure 3

SK constraints on the minimal HNL scenario at 90% CL, projected onto the plane \(|U_e|^2\) vs \(m_N\) (for \(U_{\mu ,\tau }=0\)). Our results (solid black lines) are compared to corresponding limits obtained for the T2K near detector [29], NA62 [26], E949 [25], PS191 [27, 28], and PIENU [49]. The line labeled as KENU was derived in Ref. [11] using precision measurements of leptonic decay channels of the kaon [50]

Full size image

In addition, the shaded light blue region in Fig. 2 shows the expectation for the minimal model outlined in Eq. (1), for \(m_N=250\) MeV and mixing matrix elements in the range \(|U_{e/\mu }|^2 \in [10^{-8},10^{-10}]\), \(|U_\tau |^2 \in [10^{-10}, 10^{-4}]\) (in agreement with current constraints from Ref. [9]). As readily seen from this figure, relevant constraints are expected in this case. In view of this result, next we derive constraints on the minimal scenario, assuming only one non-vanishing \(U_\alpha \) at a time. Our results are shown in Fig. 3 for \(|U_e|^2\) (left panel) and for \(|U_\mu |^2\) (right panel), at 90\(\%\) confidence level (CL). In the case of \(|U_e|^2\) the contribution from \(\pi ^\pm \) decays is clearly dominant for \(m_N < 140\) MeV, as can be seen from the peak in sensitivity at around 0.1 GeV. We find that the limits derived from our simplified analysis is already able to set tight constraints on the mixing of HNL with the electron and muon neutrino sectors, between \(10^{-6}\) and \(10^{-7}\) for \(m_N\) in the range between 150 MeV and 450 MeV. Our limits are also compared with those obtained from displaced decay searches at PS191 [27, 28] and at the T2K near detector [29], as well as from peak searches in E949 [25], PIENU [49], and NA62 [26]. We also show the resulting bound derived in Ref. [11] from the measurement of the kaon decays into electrons or muons [50]. In the case of \(|U_e|^2\), the limits obtained from SK are comparable or even better than analogous limits from peak searches, while they are not competitive with those from displaced decay searches. In the case of \(|U_\mu |^2\), peak searches in E949 also yield better constraints than our limits from SK data.

Finally, while the HNL cannot be produced via \(|U_\tau |\) in K or \(\pi \) decays (because it is not possible to produce the HNL together with a \(\tau \) lepton in this case), competitive limits can still be derived on \(U_\tau \) if we allow for non-vanishing \(|U_\alpha |,~ \alpha =e,\mu \), even if these are well below present bounds from laboratory experiments. The reason is that, in this case, the HNL could be produced via the mixing \(U_e\) or \(U_\mu \), and a large \(U_\tau \) can induce a significant decrease in its lifetime of the HNL while allowing for a significant branching ratio into e-like or \(\mu \)-like events through NC-mediated decays. Therefore, in Fig. 4 we show the sensitivity to \(|U_\tau |\), as a function of \(m_N\), for fixed \(U_e=10^{-8}\) or \(|U_\mu |=10^{-8}\) (which are both below the best present upper bounds). Our limits obtained in this way are already much better than existing direct constraints on \(|U_\tau |\) from CHARM [51], which however have been obtained assuming vanishing values for \(|U_e|\) and \(|U_\mu |\). A similar exercise can be done for PS191, which would probably lead to better limits on \(|U_\tau |\) than the SK results, for the same assumed values of \(|U_e|\) and \(|U_\mu |\).

Fig. 4
figure 4

SK constraints on the minimal HNL scenario at 90% CL, projected onto the plane \(|U_\tau |^2\) vs \(m_N\), for \(|U_{e}|^2=10^{-8}\) and \(|U_\mu |^2=10^{-8}\). Our results (black lines) are compared to the limits from CHARM [51] (which however have been obtained under the assumption \(|U_e|^2 = |U_{\mu }|^2 = 0\))

Full size image

5 Conclusions

In this letter, we have used the latest publicly available SK data to derive strong constraints on HNL production from kaon and pion decays in the atmosphere. Using a \(\chi ^2\) analysis, and binning our events only in \(\cos \theta \), we find that SK data is able to provide strong constraints on the minimal HNL scenario for masses between 150 MeV and 400 MeV. It is therefore expected that a more detailed analysis performed by the collaboration may be able to significantly improve over our results. We have also shown our limits in the Br vs \(c\tau \) plane, which is applicable to a wider range of NP models where the HNL interact with the SM not only via mixing but also through other interactions (such as, for instance, a dipole moment). Finally, we have used our results to show how, in the case of non-vanishing \(U_e\) or \(U_\mu \) well below current constraints from laboratory searches, SK data could be used to set strong bounds on \(U_\tau \), well below the direct limits presently available from CHARM data. We expect a similar improvement of PS191 bounds under similar assumptions, though.