For the simulation of the production of HNLs from light mesons, we used the data provided by the DUNE Beam Interface Working Group (BIWG) [22], which makes use of GEANT4 [23, 24] and FLUKA [25, 26]. These data include information about the decay positions and momenta of pions, kaons, and muons after they exit the focusing horns. The most abundant light parent in DUNE is the pion, followed by kaons and finally muons, as shown in Figure 4. In this work, we will consider that the neutrino CC event rates might have an overall normalization uncertainty of up to 20% due to uncertainties in the modeling of the production of mesons and leptons at the DUNE target and neutrino cross sections. We encapsulate this uncertainty by a parameter σ _a that varies from 0.05 to 0.2. Setting σ _a = 0.05 is equivalent to assume a 5% overall normalization uncertainty in the DUNE neutrino CC event rates, whereas σ _a represents an uncertainty of 20%.

The production of HNLs from heavier particles such as D mesons and τ leptons is also possible, but it is expected to have a negligible effect on the active neutrino flux, which is totally dominated by the production from lighter mesons. In order to test the relevance of HNL production from these heavy particles, we used PYTHIA8 [27] to estimate the neutrino flux generated by D 0 , D ̄ 0 , D ± , D s ± and τ ^± at DUNE. We observed that these heavy parents do not contribute significantly to the DUNE neutrino flux, and hence, the production of HNLs coming from them will have a negligible effect on the number of CC events. Consequently, our analysis is restricted only to the production of HNLs from light mesons and muons.

The production and decay chain of an HNL will depend on its mass, the mass of its parent, the nature of its parent (lepton, scalar meson, or vector meson), the parent decay channel, the HNL nature (Dirac or Majorana), the HNL decay channel, and the value of the mixing parameter involved. In principle, we could turn on, simultaneously, the three mixing parameters |U _α4|², α = e, μ, and τ; however, in our analysis, we will consider only one non-zero mixing parameter at a time.

Given the HNL mass and nature, we gave PYTHIA8 the kinematic information on the parents and let it handle the kinematics of all the HNL production and decay chain up to final active neutrinos. As expected, the HNL production and decay channels are weighted with their corresponding branching ratios.

In Figure 5, we show the number of HNLs produced at DUNE from meson decays in 1 year and in the neutrino mode for |U _μ4|² = 10^–4. The production from pion decays dominates at low masses, followed by charged and neutral kaons. The spectrum endpoint for pions and kaons corresponds to the maximum allowed m _N displayed in Table 1 when they decay into muons. For completeness, we also present the production from charmed mesons, which, as expected, is comparatively smaller and completely overshadowed for masses below 387.81 MeV. Above this threshold, HNL production from pions and kaons is kinematically forbidden, and the contribution from charmed meson decays dominates. This contribution is several orders of magnitude smaller than the one from light mesons, as we already claimed.

We focus on the active neutrinos produced from the HNL decays. We are interested in differentiating the number of these neutrinos that fall within the detector’s geometrical acceptance from those outside of it. With this aim, we parametrize the probability that an active neutrino hits the detector by two distances along the HNL propagation axis. These distances represent two different decay vertices of the HNL and are calculated considering the geometrical coverage of the detector and the kinematical information provided by PYTHIA8, which depends on its lifetime, production vertex, velocity, and the direction of the propagation of the active neutrino. The aforementioned probability is given by w d 1 , d 2 = exp − d 1 v γ τ 0 − exp − d 2 v γ τ 0 , where v is the HNL’s velocity, γ is the Lorentz factor, and τ ₀ represents proper lifetime.

For illustrative purposes, we present, in Figure 6, the scheme of the explained above, for the case when the HNL moves along the beam axis. It is clear that our analysis is general and takes into account the tridimensional shape of the LArTPC and all the possible ways in which an active neutrino might enter the detector, including cases where the HNL is outside the detector coverage.

It is important to mention that when we deactivate the HNL production, we reproduce the (pure SM) active neutrino fluxes arriving at the LArTPC predicted by the DUNE collaboration [28].

Figure 7 displays the average HNLs’ decay positions measured from the target and projected along the Z-axis for |U _μ4|² = 10^–4 and |U _μ4|² = 10^–1 and for Dirac and Majorana HNLs. The dotted line represents the position of the LArTPC, which is located at z = 574 m. Given that the lifetime of the HNL is inversely proportional to |U _μ4|², we can see that, as long as the mixing decreases, the average decay positions at Z increase. In the mass range, we studied, for |U _μ4|² = 10^–4, on average, all the HNL decays behind the LArTPC; hence, one active neutrino is lost in the DUNE flux at the LArTPC for each HNL produced. On the other hand, for |U _μ4|² = 10^–1, the average HNL decay position coincides with the LArTPC location at m _N ≈ 255 MeV, which implies that, above this mass, the HNL decay mainly before the detector.

We also note that in both cases, there is a small increase in the average decay positions at approximately 30 MeV. This happens because the production of HNLs from pion decays becomes kinematically forbidden around this energy and decays from kaons starting to dominate. This makes the average HNL more energetic, and therefore, it can travel larger distances before decaying.

The existence of HNLs forces us to modify the neutrino oscillation probabilities. Therefore, the effects of neutrino oscillations have to be taken into account in our simulations. Particularly, the place where neutrino oscillations can affect our results is in the disappearance of active neutrinos produced in meson decays. The survival probability of these active neutrinos is given by P ν α → ν α = − 4 1 − | U α 4 | 2 | U α 4 | 2 ⁡ sin 2 1.27 m N 2 L E e − Γ 4 L 2 + 2 1 − | U α 4 | 2 | U α 4 | 2 e − Γ 4 L 2 + 1 − | U α 4 | 2 2 + | U α 4 | 4 e − Γ 4 L , (6) where E represents the energy of the active neutrino, L is the distance that it travels before reaching the DUNEND, and Γ₄ is the decay rate of the HNL, and we have also considered that the mass of the active neutrino is negligible when compared to the HNL mass m _N . This survival probability will effectively decrease the number of active neutrinos that reach the DUNEND and the number of neutrino CC events at the near detector complex. For completeness, we incorporated Eq. 6 in our simulations as an extra weight for each active neutrino.

There is also the possibility of oscillation of HNLs into active neutrinos. However, since the HNL flux is very small when compared to the active neutrino flux, the effects of these oscillations in the neutrino CC event rates are negligible and were not considered in this work.

As we can infer from what we have shown before, the DUNE neutrino flux fired at the DUNEND will be affected by the production of HNLs. Each HNL produced from the decay of its parent meson (or muon) replaces one active neutrino in the SM DUNE neutrino flux. In principle, there is a possibility to recover this active neutrino since the HNL can decay into one or more active ones, which, depending on their direction, could or not impact the DUNEND. However, as demonstrated in Figure 7, it is unlikely that a relevant portion of these spurious active neutrinos would be created before or inside the LArTPC of the DUNEND for the mass range used in this work. This decrease in active neutrinos translates into a decrease in the CC event rates at the LArTPC. Our strategy is to use this deficit of CC events as an indirect signal of the existence (production) of HNLs at the DUNE neutrino flux. Hence, in that sense, we are conducting an indirect search for HNLs. This indirect method for searching HNLs is complementary to the direct searches [7], which look for HNL decays inside one of the DUNE’s detectors. As we will show in the following sections, our method can work comparatively better than direct searches for masses below 10 MeV and is sensitive to masses below 1 MeV, a region primarily inaccessible through direct searches.

The deficit in the total CC event rates depends on the mass of the HNL, the value of |U _α4|², and the off-axis position of the detector. In order to have the first estimate of the maximum significance of this deficit allowed by current limits on the mixing parameters, we calculated the active neutrino flux in the presence of HNLs using the maximum values of |U _α4|² allowed by accelerator experiments at the 90% confidence level [29] and then convoluted these fluxes with GENIE 2.8.4 [30] CC inclusive cross sections.

In order to get an idea of the significance of the change in the neutrino CC event rates, we will ignore all systematic uncertainties in the neutrino flux prediction and work in the ideal case of no systematic uncertainties σ _a = 0. Figure 8 shows the ν _μ CC event rates at the LArTPC for m _N = 1 MeV and |U _μ4|² = 10^–2 assuming Majorana neutrinos, on-axis position, and 10 years of operation (five in the neutrino mode and five in the antineutrino mode) and σ _a = 0. The significance of the change in the number of the CC events in each bin is estimated by N σ = | N BSM − N SM | N SM = | Δ N | σ , where N ^SM represents the expected number of CC events assuming only SM interactions and N ^BSM the number of CC events when HNLs are produced. As we mentioned before, we are also ignoring all normalization uncertainties in the CC event rates so that σ = N SM is the uncertainty in each bin. Due to the high luminosity of the DUNE experiment, under this setup, the production of HNLs causes a decrease in the total number of CC event rates on the order of 10⁶ events near 2.5 GeV. This implies a deviation from the SM prediction by approximately 100σ around this energy. This indicates that DUNE’s sensitivity to |U _μ4|² might be beyond the current experimental limits for this particular HNL mass.

As the HNL mass increases, its production is suppressed, and consequently, its presence on the active neutrino flux is reduced. As an example of the latter, we displayed in Figure 9 the event rates for m _N = 3 MeV and the maximum value allowed for |U _α4|² by experiments at the 90% confidence level for this mass. In this case, there is a (small) deviation, from the SM prediction, lower than 1σ. This happens because of the tighter constraint on the mixing parameter.

We have shown that in the ideal case of no systematic uncertainties σ _a = 0, DUNE will have good sensitivity for indirect hints of the existence of low mass HNLs, which are evidenced by a decrease in the neutrino CC event rates at the LArTPC. Of course, once systematic uncertainties are considered, the sensitivity and the limits are expected to decrease considerably.

In order to estimate the future sensitivity of DUNE to HNLs due to the deficit of neutrino CC events, we have to consider that the predictions of our simulations carry systematic uncertainties related to the distributions of hadron production at the DUNE target, the neutrino CC cross section uncertainties, among others. We will incorporate these uncertainties in our calculations by assuming an overall normalization uncertainty in the spectra, which, in practice, means that the values of the event rates are not completely known and can fluctuate by a certain amount. This overall normalization uncertainty will be represented by the parameter σ _a that takes the values σ _a = 0.05, 0.1, and 0.2, which are equivalent to overall normalization uncertainties of 5%, 10%, and 20%, respectively. We are also considering shape uncertainties in each bin that are represented by the parameters σ _ai ; for simplicity, we consider σ _ai = σ _a for all bins. We estimate the sensitivity of DUNE to (m _N , |U _α4|²) through the following χ ² [31]: χ 2 = a 2 σ a 2 + ∑ ν e , ν μ , ν ̄ e , ν ̄ μ ∑ i = 1 nbin a i 2 σ a i 2 + ∑ i = 1 nbin N i SM − N i BSM 1 + a + a i 2 N i SM , where N i BSM represents the neutrino CC events in the ith bin when HNLs are produced and N i SM is the DUNE prediction of CC events in the ith bin, according to the standard model. The nuisance parameters a and a _i encompass the normalization uncertainties and allow for the values of N i BSM to fluctuate; these parameters are always profiled in the calculation of χ ². We must note that the fact of combining all the neutrino flavors in our definition of χ ² is fundamental for improving the sensitivity of our results. The tau neutrinos are not considered since their contributions to χ ² are negligible.

The deficit of neutrino CC events at DUNE is an indirect signal of HNLs. Therefore, in the case that no significant deficit is found, the absence of this deficit can be used to set limits on the values of the parameters (m _N , |U _α4|²) with a particular confidence level. We calculated the value of χ ² in the parameter space 10^–6 MeV < m _N < 10⁷ MeV and 1 0 − 12 < | U α 4 | 2 < 1 for α = e, μ and σ _a = 0.5, 0.1, and 0.2 and then used these values to estimate the limits that DUNE might be able to set to the parameters (m _N , |U _α4|²) at the 90% confidence level.

Our results are presented in Figure 10. The left panel of this figure shows the estimated DUNE sensitivity to |U _μ4|² at the 90% confidence level on the LArTPC assuming Majorana neutrinos, 10 years of operation (five in neutrino and five in antineutrino modes), and on-axis position. In our analysis, the CC event rates from all neutrino flavors are considered (read the discussion at the end of Section 2). For masses close to 1 eV, the limits decrease because, for the typical energies and flight distances of active neutrinos at DUNEND, the probability of neutrino oscillations into HNLs tends to be zero as the value of m _N approaches 1 eV. Right above 1 eV, the limits start to oscillate since the survival probability of the active neutrinos is sensitive to m _N . For masses between 10 eV and 10 MeV, the limits are independent of m _N . The latter is because of three factors. The first one is the averaging out of the neutrino oscillations into HNLs for large values of m _N . The second one is that, for these very low masses, the total number of HNLs produced is practically independent of m _N (see Figure 5). The other factor is that the HNL lifetime for lower masses is enormous (see Figure 7), decaying all of them far away from the detector without the possibility of leaving a trace on it. As we already know, above m = 33.91 MeV, the production channel π ⁺ → μ ⁺ N is kinematically forbidden, and there is a sudden loss in the sensitivity. As the mass increases, production from charged kaons starts to dominate and does so up to the end of the curve, which is at 387.81 MeV. For instance, for σ _a = 0.05 and σ _a = 0.2, the sensitivity of DUNE below 10 MeV is around |U _μ4|² < 2 × 10⁻² and |U _μ4|² < 8.5 × 10⁻², respectively. We point out that even in the conservative case of σ _a = 0.2, our limits are competitive with direct searches below 1.3 MeV.

The right panel of Figure 10 shows the expected DUNE sensitivity when we turn on |U _e4|² being the other sensitivity zero. The rest of the characteristics are the same as for the left panel. In general, the sensitivity pattern is similar to the pattern observed for the left panel. The limits oscillate close to 1 eV, and for higher masses, they become mass-independent since most HNLs decay behind the LArTPC. Above 10 MeV, the pion decay channel π ^± → e ^± N starts to dominate because, in contrast to π ± → e ± ν e ( − ) , it is less suppressed by helicity due to the larger size of the HNL mass. This effect decreases the number of both ν _e and ν _μ CC events according to the suppression factor shown in Eq. 5, affecting the CC event rates of both electron and muon neutrinos. At approximately 139 MeV, the HNL production from pion decays becomes kinematically forbidden, which translates into a decrease in the sensitivity. Finally, the curve ends when the production from kaons is kinematically forbidden at 493.17 MeV. For σ _a = 0.05 and σ _a = 0.2, the sensitivity of DUNE below 10 MeV is around |U _e4|² < 1.5 × 10⁻² and |U _e4|² < 6.5 × 10⁻², respectively. Even in the conservative case of σ _a = 0.2, our limits are competitive with direct searches below 1.3 MeV and also provide a small increase in sensitivity by a factor of 1.5 at approximately 5 eV in comparison with experimental constraints.

Although we are making our calculations for 10 years of exposure, it is important to point out that our sensitivity for |U|² increases only slightly when compared with 1 year of exposure. If we had not included systematic uncertainties, the limits would roughly improve as | U | 2 T , where T represents the exposure time; in this ideal scenario, after 10 years of operation, the limits would improve by a factor of around 1 / 10 ≈ 0.32 . However, introducing uncertainties in our χ ² prescription heavily penalizes the sensitivity of our approach: in this more realistic scenario, after 10 years of operation, the limits improve by only a factor of approximately 0.9 in comparison with 1 year of exposure. Therefore, in the context of our analysis, the first year of operation of DUNE is the most important.

Another important remark must be done about the effects of neutrino oscillations in this work. Neutrino oscillations involving HNLs are only relevant when m _N ∼1 eV. Since our analysis starts approximately at 1 eV, the effects of neutrino oscillations will only be visible as a wiggle at the beginning of our sensitivity plots. For completeness, in Figure 11, we show a zoom of the left plot in Figure 10. We can see that close to 1 eV, the sensitivity oscillates as expected, but this effect is small and only restricted to the low tail of our sensitivity plot.

We must point out that our results are blind to the Dirac or Majorana nature of the HNL. The distinction between Dirac and Majorana HNLs is usually performed in direct searches by analyzing the distributions of charged mesons and leptons produced when the HNL decays inside the detector. We are not looking into the direct search mode since it has already been discussed in [7]. In addition to their decay products, Dirac and Majorana HNLs can also be differentiated by their lifetimes due to the factor of two present in Eq. 5. However, this effect is not relevant for us because, for the mass range we studied and small mixings, almost all the HNL decays occur behind the LArTPC, as shown in Figure 7. Furthermore, as we have discussed in Section 2, for very low m _N , the Dirac and Majorana neutrinos are indistinguishable. Thus, we can conclude that nearly all the active neutrinos produced from the HNL decays are lost independently of the nature of neutrinos. In this way, the critical magnitude in our analysis is the production rate of HNLs, which is independent of the nature of neutrinos, so the deficit of the CC event rates is independent too. Therefore, it would not be possible to distinguish between Dirac and Majorana neutrinos through the approach presented here.

The DUNE experiment also considers the possibility of moving the DUNE near detectors horizontally, a setup known as DUNE-PRISM. We move the LArTPC by up to 30 m horizontally while maintaining the rest of the simulation parameters and study the impact in our estimated sensitivities. The results are shown in Figure 12, where all the lines represent the sensitivities at the 90% confidence level and the dashed curves represent the on-axis sensitivities. We see that the effect of moving the detector to an off-axis position does not affect considerably the limits, although the curves are less smooth due to the decrease in statistics. However, we see that the sensitivities increase at off-axis positions for masses close to 100 MeV.

We also explore the potential to constrain the (m _N , |U _α4|²) parameter space region in the context of this indirect search. So, assuming that the disappearance CC events are originated by the presence of HNLs within the neutrino beam, we perform χ ² analysis fixing our simulation in certain values of (m _N , |U _α4|²). The 95% confidence regions for m _N = 0.1 MeV and |U _μ4|² = 5 × 10⁻² are presented in Figure 13 for σ _a = 0.05 (red), σ _a = 0.1 (blue), and σ _a = 0.2 (green). For the case of 20% uncertainty, the analysis only gives upper limits to the values of |U _μ4|². For 5% and 10% uncertainties, the regions are bounded to the right but extend to the left up to m _N = 1 eV, a mass degeneracy that reflects the fact that our approach is not sensitive to m _N for low masses. For the case σ _a = 0.05 (red), the 95% confidence region is sufficiently small that it is possible to constraint |U _μ4|² within an uncertainty of 50%. However, when we include larger systematic uncertainties such as σ _a = 0.1 (blue), we found that we can only constrain the value of |U _μ4|² within an uncertainty of 100%.