New constraints on heavy neutral leptons from Super-Kamiokande data

Heavy neutral leptons are predicted in many extensions of the Standard Model with massive neutrinos. If kinematically accessible, they can be copiously produced from kaon and pion decays in atmospheric showers, and subsequently decay inside large neutrino detectors. We perform a search for these long-lived particles using Super-Kamiokande multi-GeV neutrino data and derive stringent limits on the mixing with electron, muon and tau neutrinos as a function of the long-lived particle mass. We also present the limits on the branching ratio versus lifetime plane, which are helpful in determining the constraints in non-minimal models where the heavy neutral leptons have new interactions with the Standard Model.


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 : a e-mail: pilar.coloma@ific.uv.es (corresponding author) where˜ ≡ iσ 2 * is the complex conjugate of the Higgs field , L α is the SM lepton doublet with flavor α, 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 ν ∝ m −1 N (Y v) 2 , and two or more heavy neutral leptons (HNL) with masses ∝ 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 α 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 α | ≡ U α 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 (τ ) depends strongly on its mass and on its mixing with the light states, in the minimal model described above it ranges between cτ ∼ (10 −4 −50)× |U α | −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 Ice-Cube [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 α | versus m N ), but also on the plane Br(K /π → N ) versus cτ , 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].

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 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 ) can be obtained as: where dn/d E stands for the differential distribution of the HNL energies while P is the parent meson flux and, in this work, P = K ± , π ± . 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 the 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), β P and γ P are the parent boost factors while τ P is its lifetime in the rest frame. For the dominant two-body decay K ± , π ± → Nl ± α (denoted generally as P → N Y ), the differential distribution dn/d E reads where Br stands for branching ratio, and Finally, the kinematical limits for E P m P are: .
(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 τ 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.
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 νN → N X, where N stands for a nucleon and X is a hadronic shower). This contribution can be estimated as follows. Assuming a flux of standard neutrinos φ ν , the HNL produced at a distance r from SK is: where σ ν→N is the HNL production cross section (which can be estimated as the standard ν NC scattering cross section multiplied by the corresponding mixing |U α | 2 ), while N A is the Avogadro number and ρ ⊕ is the Earth density at distance 1 While significant variations are expected for the prompt neutrino flux if the cosmic ray or hadronic interaction models are changed (see e.g. Ref. [43]), the conventional neutrino contribution is understood at the O(10 − 20%) level. We expect similar variations in our results, if these assumptions are modified. In this work, we consider only π and K decays. For comparison we also show the corresponding result for twobody decays of heavier mesons (D, D s ) and leptonic (three-body) τ 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 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.

HNL decays in Super-Kamiokande
The flux that arrives to the detector for HNL produced in meson decays, decay N , is obtained integrating over all LLPs produced at different distances , weighted by their corresponding survival probabilities, as where L decay is the decay length of the HNL in the laboratory frame. The maximum distance max ≡ (h max , θ) is a simple function of the maximum height of the atmosphere where cosmic showers are produced, h max 80 km, and the zenith angle. Analogously, the flux entering the detector for HNL produced in the interaction of SM neutrinos in the Earth, int N is: where R max (θ ) = 2R ⊕ cos θ is the maximum distance traveled through the Earth for trajectories with a zenith angle θ , R ⊕ being the Earth's radius.
where N ≡ int N + dec N , and A 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: Here 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τ but, within the range cτ ∈ [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 μand e-like fullycontained multi-GeV events, adding the single-and multiring samples together (labeled as "multi-GeV" and "multiring" in Fig. 5 in Ref. [48], respectively). In the case of e-like events, we add both the ν e andν 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 θ . 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 vari-ables. 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 α = 0.75 both for μ and e-like, in line with the values quoted in Ref. [48] for the multi-GeV ν 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 θ 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 θ in a given sample can therefore be computed as: where cos θ min i and cos θ max i are the lower and upper limits of the bin. Here, Br(α-like) stands for the total branching ratio for all decay channels including muons, electrons or photons in the final state, depending on the sample (α) considered. For example, in the case of μ-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 → νπ 0 , since the π 0 decays promptly to two photons).

Results
In this section we derive limits on HNL production by performing a χ 2 fit to the SK data. A Poissonian χ 2 function has been used: where the sum runs over the angular bins, and α = {elike, μ-like}. Here, n α i stands for the data observed in each bin while N α i is the predicted number of signal events and B α i is the background prediction, which includes the predicted number of atmospheric neutrino events in the SM. In our χ 2 calculations, we consider separately the e-and μ-like samples and we add their two contributions to the total χ 2 . Therefore, our limits will be derived using 20 degrees of freedom, corresponding to the total number of bins in cos θ . In our calculations we find, however, that the sensitivity is largely dominated by the e-like contribution since the size of the branching ratio Br(e-like) is much larger than Br(μ-like) in the mass range considered.
First we show in Fig. 2 the results on the plane Br(K /π → N ) × Br(N → visible), where Br(N → visible) accounts for the probability that the HNL decays visibly in the detector. Our results are presented as a function of cτ , 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 → νγ 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 ± decay (pink lines) fall short to probe this region. However, SK would be sensitive to non-minimal models with HNL produced in π ± 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.
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/μ | 2 ∈ [10 −8 , 10 −10 ], |U τ | 2 ∈ [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 α at a time. Our results are shown in Fig. 3 for |U e | 2 (left panel) and for |U μ | 2 (right panel), at 90% confidence level (CL). In the case of |U e | 2 the contribution from π ± 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 μ | 2 , peak searches in E949 also yield better constraints than our limits from SK data. Finally, while the HNL cannot be produced via |U τ | in K or π decays (because it is not possible to produce the HNL together with a τ lepton in this case), competitive limits can still be derived on U τ if we allow for non-vanishing |U α |, α = e, μ, 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 μ , and a large U τ can induce a significant decrease in its lifetime of the HNL while allowing for a significant branching ratio into e-like or μ-like events through NC-mediated decays. Therefore, in Fig. 4 we show the sensitivity to |U τ |, as a function of m N , for fixed U e = 10 −8 or |U μ | = 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 τ | from CHARM [51], which however have been obtained assuming vanishing values for |U e | and |U μ |. A similar exercise can be done for PS191, which would probably lead to better limits on |U τ | than the SK results, for the same assumed values of |U e | and |U μ |.

Conclusions
In this letter, we have used the latest publicly available SK data to derive strong constraints on HNL production from  [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]  Our results (black lines) are compared to the limits from CHARM [51] (which however have been obtained under the assumption |U e | 2 = |U μ | 2 = 0) kaon and pion decays in the atmosphere. Using a χ 2 analysis, and binning our events only in cos θ , 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τ 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 μ well below current constraints from laboratory searches, SK data could be used to set strong bounds on U τ , well below the direct limits presently available from CHARM data. We expect a similar improvement of PS191 bounds under similar assumptions, though. MSCA-RISE-2015. The work of I.M.S. is supported by the U.S. Department of Energy under the award number DE-SC0020250.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: Data sharing is not applicable to this article, as no datasets were generated or analysed during this study.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .