Indirect search of heavy neutral leptons using the DUNE near detector

We evaluate the potential of the DUNE near detector (DUNEND) for establishing bounds for heavy neutral leptons (HNLs). This is achieved by studying how the presence of HNLs affects the production rates of active neutrinos, therefore creating a deficit in the neutrino charged current (CC) events at the LArTPC of the DUNEND. The estimated bounds on HNLs are calculated for masses between 1 eV and 500 MeV. We consider 10 years of operation (five in neutrino and antineutrino modes) and obtain limits of |U_μ4|² < 9 × 10⁻³ (4 × 10⁻²) and |U_e4|² < 7 × 10⁻³ (3 × 10⁻²) for masses below 10 MeV and a 5% (20%) overall normalization uncertainty in the neutrino charged current event rate prediction. These limits, within the region of masses below 2 (10) MeV, are better than those that can be achieved by DUNE direct searches for the case of a 5% (20%) uncertainty. When a conservative 20% uncertainty is present, our limits can only improve current constraints on |U_e4|² by up to a factor of 3 in a small region of approximately 5 eV and set limits on |U_μ4|² in a mass region free of constraints (40 eV–1 MeV).

1 Introduction

Heavy neutral leptons (HNLs) are singlet (right-handed) fermion states introduced for explaining the non-zero neutrino masses, interacting via Yukawa coupling with the Higgs boson and the leptonic doublet, a Dirac mass term, and also appearing in the Majorana mass term. The nearly sterile states that arise after the diagonalization of the mass terms interact with matter via suppressed mixing to the active neutrinos of the Standard Model (SM) [1, 2].

The HNLs are candidates to solve important particle physics and cosmological issues [1]. They can help explain the smallness of the active neutrino masses via the seesaw mechanism [3], act as possible dark matter candidates [4], and also explain the baryon asymmetry of the universe through their role in leptogenesis (see [1, 5] and references therein). On the other hand, neutrino oscillations involving light sterile states have been proposed to explain the excess of electron antineutrino and neutrino events at LSND and MiniBoonE, respectively, as well as the deficit of electron antineutrino events at reactor experiments [6]. The HNL masses required for solving the previously mentioned problems fall within a mass range that spans from keV to TeV. As a consequence of their relevance, there have been several HNL searches in this wide mass range, placing limits on the possible values of the HNL mass m_N and its mixing to the SM neutrinos |U_α4|² [2, 7].

In particular, searches of HNLs in the range of masses of 1–400 MeV have been conducted in accelerator-based experiments through searches for low-energy peaks in the energy spectrum of the muons resulting from pion (π^± → μ^±ν_H) and kaon decays (K^± → μ^±ν_H) [8–11]. With no positive results found so far, they obtain upper bounds for |U_μ4|² such as 10^–6 for m_N ∼ 10 MeV and 10^–9 for m_N ∼ 300 MeV.

This work aims to assess the sensitivity of the DUNE in setting upper limits for |U_μ4|² and |U_e4|² for masses below 500 MeV. We achieved the latter by comparing the energy distributions of the neutrino CC event rates with and without HNLs at the DUNE near detector (DUNEND) [12]. We found that the presence of HNLs creates a deficit of CC events that is not generated by neutrino oscillations but instead by the combination of kinematic effects in the production and decay chains of HNLs: the decrease in the branching ratios of active neutrino production, the large lifetimes of HNLs, and the fact that active neutrinos born from HNLs have angular distributions spanned outside the detector coverage. We consider the decrease in CC events as an indirect signal of HNLs and use it to set limits on the mixing parameters. Additionally, we present an analysis of the possibility of finding confidence regions for the values of (m_N, |U_α4|²) if a deficit of CC events is found at DUNE [13].

This paper is summarized as follows: in Section 2, we discuss the theoretical framework of HNL production and decay. Then, in Section 3, we describe the experimental setup. In Section 4, the details of our simulation are given, while in Section 5, our results are presented. We draw our conclusion in the final section.

2 Theoretical framework

As we already mentioned, the nearly sterile mass eigenstates couple to the active flavor states via an extended version of the Pontecorvo–Maki–Nakagawa matrix (PMNS) [14], which can be expressed as follows:

${\nu }_{\alpha }={\displaystyle \sum _{i=\mathrm{1,2,3}}}{U}_{\alpha i}{\nu }_i+U_{\alpha 4}N,$

where N represents the HNL field. It is also helpful to write the new active neutrino flavor states in terms of the flavor states of the SM ${\nu }_{\alpha }^{\text{SM}}$, which represent the neutrino flavor states when the values of the 3 × 3 PMNS mixing matrix are assumed. This can be done by the approximation [15].

${\nu }_{\alpha }\approx {\nu }_{\alpha }^{\text{SM}}\left(1-\frac{|U_{\alpha 4}{|}^2}{2}\right)+U_{\alpha 4}N.$

Due to the connection mentioned above, the HNLs can be produced in any weak decay involving active neutrinos. The production rate of HNLs depends kinematically on its mass m_N, the strength of its mixing to active neutrinos |U_α4|², and the nature of the decaying particle that produces it, which, from now on, is referred to as its parent. In this work, we are interested in HNLs with masses below the kaon mass (m_K). The production of HNLs from kaon and pion decays, followed by the muon decays, dominates at the typical energies of beam dump experiments such as DUNE. Their production from heavier particles, such as D mesons or τ leptons, is also possible, but it is rare since the production of the latter is heavily suppressed in comparison to the light mesons. Table 1 shows the dominant HNL production channels from light leptons and mesons, along with the maximum kinematically allowed values of the masses for the HNLs. A rough estimation of these values is obtained by subtracting the total rest mass of the particles produced, other than the HNLs, from the corresponding mass of their parent particles.

Table 1

Table 1. Channels considered for the production of HNLs. The maximum possible value of m_N is shown for each channel. Charged conjugate channels were also considered.

We calculated the branching ratios for HNL production by using the formulas from [16]. For instance, Figure 1 shows the branching ratios of the dominant HNL production channels below the kaon mass for |U_μ4|² = 1. We can note that almost all the branching ratios decrease with m_N, with the only exception being the leptonic decays of charged kaons, K^± → Nμ^±. Above 34 MeV, the production from pions is kinematically forbidden; this is important since this means that all heavy neutral leptons above this mass will be produced only from kaon decays. As the value of m_N increases, the branching ratio of K^± → Nμ^± keeps increasing as well, surpassing the branching ratios of K^± → Nπ⁰μ^± at approximately 80 MeV and $K_L^0\to N{\pi }^0{\mu }^{\pm }$ at approximately 160 MeV. Finally, the branching ratio of K^± → Nμ^± reaches its maximum at approximately 260 MeV and then decreases until it is kinematically forbidden. The endpoint of each branching ratio corresponds to the maximum m_N given in Table 1.

Figure 1

Figure 1. Branching ratios of the dominant HNL production channels for |U_μ4|² = 1.

The production of HNLs via semileptonic decays involves hadronic currents that cannot be calculated from first principles due to the non-perturbative nature of QCD at low energies. Therefore, the dynamics of these decays are modeled by form factors that represent the momentum distribution of the quarks inside the mesons and parametrize the momentum transfer between the hadronic current and the lepton pair [17]. For all the semileptonic decays shown in Table 1, we used the form factors presented in [16].

After their production, all the HNLs propagate and then decay on flight via mixing with active neutrinos. Table 2 shows all the decay channels for the HNLs considered in this work. We included all the kinematically allowed decays to final states involving pseudoscalar mesons and pure leptonic decays for m_N < m_K. A more complete table is found in [18].

Table 2

Table 2. HNL decay channels considered in this work. The minimum required value of m_N is shown for each channel.

The partial width of an HNL decay channel involving a final lepton l_α or light neutrino ν_α is directly proportional to the mixing parameter squared |U_α4|². Therefore, the total width and lifetime of the HNLs also depend on the relevant mixing parameters. The lifetime dependence on the values of |U_α4|² can have a huge impact on the position of the decay vertex of the HNL and, hence, on its possible signal in the detector. Setting small values for the |U_α4|² means that HNLs are being produced at a lower rate, but at the same time, these HNLs have a greater lifetime and therefore decay further away from the detector.

When we determine the individual partial widths of each channel, there is a factor of two that differentiates between the decays of Dirac and Majorana HNLs [18]. For instance, a Dirac HNL can decay to charged pions only via N → e⁻π⁺, while a Majorana HNL can also decay through N → e⁺π⁻. This evidently has an effect on the rates of π⁺/π⁻ production from HNL decays but does not affect the partial decay widths. This means that CC-mediated channels have the same partial widths for Dirac and Majorana neutrinos:

$\begin{array}{c}\hfill \mathrm{\Gamma }\left(N_M\to l^{-}{X}^{+}\right)=\mathrm{\Gamma }\left(N_D\to l^{-}{X}^{+}\right),\\ \hfill \mathrm{\Gamma }\left(N_M\to l^{+}{X}^{-}\right)=\mathrm{\Gamma }\left({\bar{N}}_D\to l^{+}{X}^{-}\right).\end{array}(1)$

On the other hand, NC-mediated channels do distinguish between Dirac and Majorana HNLs. This is because the contractions of the NC operator add an additional contribution to differential decay widths of the Majorana HNLs [18, 19].

$d\mathrm{\Gamma }\left(N_M\to \nu X\right)=d\mathrm{\Gamma }\left(N_D\to \nu X\right)+d\mathrm{\Gamma }\left({\bar{N}}_D\to \bar{\nu }X\right).$

Therefore, a factor of two appears when comparing the partial widths of NC-mediated decays.

$\mathrm{\Gamma }\left(N_M\to \nu X\right)=2\mathrm{\Gamma }\left(N_D\to \nu X\right).(2)$

Equations 1 and 2 imply that the total widths (Γ_T) of Majorana and Dirac HNLs are related by

${\mathrm{\Gamma }}_{\text{T}}\left(N_M\right)=2{\mathrm{\Gamma }}_T\left(N_D\right),$

which translates into a difference between their lifetimes:

$\tau \left(N_M\right)=\frac{1}{2}\tau \left(N_D\right).$

For very low masses (m_N ≪ m_e), the factor of two in Eq. 2 disappears [20], making the total widths and lifetimes of Dirac and Majorana HNLs indistinguishable. Part of the mass range that we will explore in this work falls in the region of very low masses.

At the end of this section, we will describe how the active neutrino flux is affected by the production of HNLs. For this purpose, we will show how the SM parent meson’s branching ratios are modified when the production of HNL occurs. Let us start by defining the SM total decay rate of the pion $({\mathrm{\Gamma }}_{\pi }^{\text{SM}})$:

${\mathrm{\Gamma }}_{\pi }^{\text{SM}}={\mathrm{\Gamma }}^{\text{SM}}\left(\pi \to e{\nu }_e\right)+{\mathrm{\Gamma }}^{\text{SM}}\left(\pi \to \mu {\nu }_{\mu }\right),$

and the decay rate with heavy neutral leptons $({\mathrm{\Gamma }}_{\pi }^{\text{BSM}})$:

$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_{\pi }^{\text{BSM}}=& {\mathrm{\Gamma }}^{\text{BSM}}\left(\pi \to e{\nu }_e\right)+{\mathrm{\Gamma }}^{\text{BSM}}\left(\pi \to \mu {\nu }_{\mu }\right)+\mathrm{\Gamma }\left(\pi \to NX\right)\hfill \\ \hfill \approx & {\mathrm{\Gamma }}^{\text{SM}}\left(\pi \to e{\nu }_e\right)\left(1-\frac{|U_{e4}{|}^2}{2}\right)\hfill \\ \hfill & +{\mathrm{\Gamma }}^{\text{SM}}\left(\pi \to \mu {\nu }_{\mu }\right)\left(1-\frac{|U_{\mu 4}{|}^2}{2}\right)\hfill \\ \hfill & +\mathrm{\Gamma }\left(\pi \to NX\right).\hfill \end{array}$

The branching ratio of ν_μ production from pion decays in the presence of HNLs can then be written as

$\begin{array}{cc}\hfill {\text{BR}}^{\text{BSM}}\left(\pi \to \mu {\nu }_{\mu }\right)& =\frac{{\mathrm{\Gamma }}^{\text{BSM}}\left(\pi \to \mu {\nu }_{\mu }\right)}{{\mathrm{\Gamma }}_{\pi }^{\text{BSM}}}\hfill \\ \hfill & \approx \frac{{\mathrm{\Gamma }}^{\text{SM}}\left(\pi \to \mu {\nu }_{\mu }\right)\left(1-\frac{|U_{\mu 4}{|}^2}{2}\right)}{{\mathrm{\Gamma }}_{\pi }^{\text{SM}}}\cdot \frac{{\mathrm{\Gamma }}_{\pi }^{\text{SM}}}{{\mathrm{\Gamma }}_{\pi }^{\text{BSM}}}\hfill \\ \hfill & \approx {\text{BR}}^{\text{SM}}\left(\pi \to \mu {\nu }_{\mu }\right)\cdot \frac{{\mathrm{\Gamma }}_{\pi }^{\text{SM}}}{{\mathrm{\Gamma }}_{\pi }^{\text{BSM}}}\left(1-\frac{|U_{\mu 4}{|}^2}{2}\right).\hfill \end{array}(3)$

A similar relation can be found for the branching ratio of ν_e production from pion decays:

$\begin{array}{c}\hfill {\text{BR}}^{\text{BSM}}\left(\pi \to e{\nu }_e\right)\approx {\text{BR}}^{\text{SM}}\left(\pi \to e{\nu }_e\right)\cdot \frac{{\mathrm{\Gamma }}_{\pi }^{\text{SM}}}{{\mathrm{\Gamma }}_{\pi }^{\text{BSM}}}\left(1-\frac{|U_{e4}{|}^2}{2}\right),\end{array}(4)$

where BR^SM(π → μ(e)ν_μ(e)) represents the branching ratio of ν_μ(ν_e) production from pion decays in the SM. We can see that the introduction of HNLs causes the production of either muon or electron neutrinos from pions to be suppressed by the factor:

${\mathcal{K}}_{\pi }^{\alpha }\left(m_N,|U_{\alpha 4}{|}^2\right)=\frac{{\mathrm{\Gamma }}_{\pi }^{\text{SM}}}{{\mathrm{\Gamma }}_{\pi }^{\text{BSM}}}\left(1-\frac{|U_{\alpha 4}{|}^2}{2}\right),(5)$

with α = e, μ. Figure 2 illustrates the dependence on m_N of the factor ${\mathcal{K}}^{\mu }$ for several parents assuming |U_μ4|² = 10^–4. For each meson, the suppression factor acts only up to a maximum HNL mass due to kinematical constraints, which are the same constraints shown in Table 1; Figure 1. Although the effect is small, the high luminosity of DUNE makes it possible to use this effect to set limits on the heavy neutral lepton parameters.

Figure 2

Figure 2. Suppression factor ${\mathcal{K}}^{\mu }\left(m_N,|U_{\mu 4}{|}^2=10^{-4}\right)$ of muon neutrino production as a function of m_N.

Thus, each particle capable of producing active neutrinos can now produce HNLs, leading to a suppression of active neutrinos at DUNE. The latter happens for all neutrino flavors even when only one mixing |U_α4|² is turned on. In fact, we can see from Eqs 3, 4 that, if we set |U_α4|² = 0, the production of the active neutrinos ν_α is still suppressed by the factor ${\mathrm{\Gamma }}_{\pi }^{\text{SM}}/{\mathrm{\Gamma }}_{\pi }^{\text{BSM}}$. As we will show further ahead, the reduction in the active neutrino flux would imply the possibility that they do not reach the DUNEND, decreasing the number of expected CC events at this facility.

3 Experimental setup

In order to simulate how the presence of HNLs affects the number of ν CC events at DUNE, we based our experimental setup in the DUNE near detector, as described in [13].

We assume that the LBNF-DUNE beam collides protons with approximately 120 GeV of energy into a graphite target, producing 1.47 × 10²¹ POTs per year. At each collision, several mesons are produced, including mostly pions, kaons, and charmed mesons.

The muons and long-lived charged mesons (π^± and K^±) produced are deflected by focusing magnetic horns located right after the target; as a consequence, their trajectories end up preferably oriented along the beam axis, as shown schematically in Figure 3. On the other hand, the trajectories of neutral mesons (D⁰, $K_L^0$ and π⁰), tau leptons, and short-lived charged heavy mesons (D^± and $D_s^{\pm }$) are not affected by the focusing horns. Most particles decay in flight inside the decay pipe, a cylinder with a length of 230 m, and a diameter of 2 m; however, a small number of long-lived particles reach the end of the decay pipe and decay at rest at the decay pipe’s surface.

Figure 3

Figure 3. Experimental setup for the LArTPC in the neutrino mode (not to scale). Charged particles are deflected by the magnetic horns.

The near detector liquid argon time projection chamber (LArTPC) is located at 574 m from the target. It is parallelepiped with a width and height (both transverse to the beam direction) of 7 m and 3 m, respectively, and a length of 5 m in the beam direction. The LArTPC is filled with a fiducial mass of 50 tons of liquid argon. There is also the multi-purpose detector (MPD), which is a magnetic spectrometer designed to study particles exiting the LArTPC and contains a one-ton high-pressure cylindrical gaseous argon time projection chamber. Since we are interested in the effects of HNLs on the ν CC events at the DUNE near detector, we will not take into account the MPD in our simulation setup because its impact on our results is negligible.

We also take into account the possibility of moving the detectors to several off-axis positions along the x-axis, a setup known as DUNE-PRISM [21].

4 Simulation route for HNLs

4.1 Parent production

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%.

Figure 4

Figure 4. Spectra of light particles capable of producing HNLs in the DUNE beam. Different bin widths have been used for different particles.

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,{\bar{D}}^0,D^{\pm },D_s^{\pm }$ 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.

4.2 Production of HNLs

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.

Figure 5

Figure 5. Heavy neutral leptons produced from mesons in 1 year in the neutrino mode for |U_μ4|² = 10^–4.

4.3 Decay of HNLs—active neutrinos

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\left(d_1,d_2\right)=\exp \left(-\frac{d_1}{v\gamma {\tau }_0}\right)-\exp \left(-\frac{d_2}{v\gamma {\tau }_0}\right),$

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.

Figure 6

Figure 6. HNL N propagates and decays into an active neutrino ν_μ. If the HNL decays between positions 1 and 2, the active neutrino ν_μ hits the LArTPC.

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.

Figure 7

Figure 7. Average HNL’s decay positions protected along the Z-axis for |U_μ4|² = 10^–4 and |U_μ4|² = 10^–1. The dotted line represents the position of LArTPC.

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.

4.4 Oscillation effects in active neutrinos from meson decays

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

$\begin{array}{ll}\hfill P_{{\nu }_{\alpha }\to {\nu }_{\alpha }}& =-4\left(1-|U_{\alpha 4}{|}^2\right)|U_{\alpha 4}{|}^2{\sin }^2\left(\frac{1.27m_N^{2}L}{E}\right)e^{-\frac{{\mathrm{\Gamma }}_{4}L}{2}}\hfill \\ \hfill & +2\left(1-|U_{\alpha 4}{|}^2\right)|U_{\alpha 4}{|}^2{e}^{-\frac{{\mathrm{\Gamma }}_{4}L}{2}}\hfill \\ \hfill & +{\left(1-|U_{\alpha 4}{|}^2\right)}^2+|U_{\alpha 4}{|}^4{e}^{-{\mathrm{\Gamma }}_{4}L},\hfill \end{array}(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.

5 Results

5.1 Impact on CC events at DUNEND

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_{\sigma }=\frac{|N^{\text{BSM}}-N^{\text{SM}}|}{\sqrt{N^{\text{SM}}}}=\frac{|\mathrm{\Delta }N|}{\sigma },$

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 $\sigma =\sqrt{N^{\text{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.

Figure 8

Figure 8. ν_μ CC event rates for m_N = 1 MeV assuming the maximum value allowed for |U_μ4|² at the 90% confidence level, on-axis position, 10 years of operation, and σ_a = 0. The error bars are amplified by 100.

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.

Figure 9

Figure 9. ν_μ CC event rates for m_N = 3 MeV assuming the maximum value allowed for |U_μ4|² at the 90% confidence level, on-axis position, 10 years of operation, and σ_a = 0.

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.

5.2 Sensitivity

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]:

${\chi }^2=\frac{a^2}{{\sigma }_a^2}+{\displaystyle \sum _{{\nu }_e,{\nu }_{\mu },{\bar{\nu }}_e,{\bar{\nu }}_{\mu }}}\left[{\displaystyle \sum _{i=1}^{\text{nbin}}}\frac{a_i^2}{{\sigma }_{ai}^2}+{\displaystyle \sum _{i=1}^{\text{nbin}}}\frac{{\left(N_i^{\text{SM}}-N_i^{\text{BSM}}\left(1+a+a_i\right)\right)}^2}{N_i^{\text{SM}}}\right],$

where $N_i^{\text{BSM}}$ represents the neutrino CC events in the ith bin when HNLs are produced and $N_i^{\text{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^{\text{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 $10^{-12}<|U_{{\alpha }^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.

Figure 10

Figure 10. Estimated limits of DUNE to |U_μ4|² (left, red) and |U_e4|² (right, blue) at the 90% confidence level by CC event disappearance at the LArTPC of the DUNEND, for 10 years of operation (five in neutrino and five in antineutrino modes) and the on-axis position. The regions of experimental constraints (gray) were taken from [29, 32, 33]. The estimated sensitivity of DUNE obtained in [7] by direct searches of HNL decays is shown for comparison.

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 ${\pi }^{\pm }\to e^{\pm }\stackrel{(-)}{{\nu }_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 $\frac{|U{|}^2}{\sqrt{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/\sqrt{10}\approx 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.

Figure 11

Figure 11. Zoom of the sensitivity to |U_μ4|² of Figure 10. The oscillation of the sensitivity near 1 eV is produced by the oscillation effects of Eq. 6.

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.

5.3 Off-axis sensitivity

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.

Figure 12

Figure 12. Comparison between on-axis (dashed) and 30 m off-axis (solid) estimated sensitivities of DUNE to |U_μ4|² at the 90% confidence level by neutrino CC event disappearance for 10 years of operation (five in neutrino and five in antineutrino modes). The regions of experimental constraints were taken from [29, 32, 33].

5.4 Allowed regions for (m_N, |U_α4|²)

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%.

Figure 13

Figure 13. Regions having 95% confidence level for m_N = 0.1 MeV, |U_μ4|² = 5 × 10⁻², 10 years of operation (five in neutrino and five in antineutrino modes), on-axis position, and several values of σ_a.

6 Conclusion

The cornerstone of this work is the analysis of other ways in which the active neutrino flux at DUNE is affected by the production of HNLs besides neutrino oscillations. We found that, in the presence of HLNs, the production rates of active neutrinos decrease and their angular distributions widen, which translates into a decrease in the number of ν_μ and ν_e CC events in the LArTPC of the DUNEND. This neutrino disappearance represents an indirect signal of HNLs at DUNE that is not due to neutrino oscillations but rather to the kinematics of the meson and HNL decays. When combined with the effects of neutrino oscillations, it is possible to use this deficit in CC event rates to estimate limits of |U_e4|² and |U_μ4|². We found that these limits are very sensitive to the uncertainty of the neutrino flux prediction at the DUNEND. In order to get conservative estimates of these limits, we considered overall normalization uncertainties of up to 20%.

For 5 year per mode (neutrino/antineutrino), on-axis configuration, and a 5% overall normalization uncertainty, we obtain limits of |U_μ4|² < 2 × 10⁻² and |U_e4|² < 1.5 × 10⁻² below 1.5 MeV. We also included a more pessimistic scenario of a 20% systematic uncertainty and were still able to set bounds of |U_μ4|² < 8.5 × 10⁻² and |U_e4|² < 6.5 × 10⁻² below 1.4 MeV. These limits are better than the limits predicted by DUNE direct searches or even placed in mass regions inaccessible to them. These bounds are still competitive for the off-axis configuration. In addition, we explore the capacity of determining the allowed parameter space region (m_N, |U_α4|²) for the specific parameter values m_N = 0.1 MeV and |U_μ4|² = 5 × 10⁻² and found that although there is a large degeneracy in the value of m_N, it is possible to constrain |U_μ4|² with uncertainties in the order of 50 (100)% for a 5 (10)% overall normalization uncertainty in the CC event rates. Finally, it is worth noting that the disappearance of CC events as an HNL signature is complementary to the direct observation or HNL decays, showing an attractive potential to be used in neutrino ear detectors with high ν CC event rates.

Data availability statement

Publicly available datasets were analyzed in this study. These data can be found at: https://glaucus.crc.nd.edu/DUNEFluxes/.

Author contributions

SC: conceptualization, data curation, formal analysis, funding acquisition, investigation, methodology, project administration, resources, software, supervision, validation, visualization, writing–original draft, and writing–review and editing. AG: conceptualization, data curation, formal analysis, funding acquisition, investigation, methodology, project administration, resources, software, supervision, validation, visualization, writing–original draft, and writing–review and editing.

Funding

The authors declare that financial support was received for the research, authorship, and/or publication of this article. AG acknowledges funding by the Dirección de Gestión de la Investigación at PUCP, through grant nos DGI-2017-3-0019 and DGI 2019-3-0044. SC acknowledges CONCYTEC for the graduate fellowship under Grant No. 236-2015-FONDECYT.

Acknowledgments

The authors also want to thank Zarko Pavlovic for useful remarks regarding the interpretation of the BIWG data, Margot Delgado de la Flor for her help at the initial stages of this work, and R. E. Shrock, C. Argüelles, I. Shoemaker, A. de Roeck, O. Peres, and W. Rodejohann for their useful remarks on our first manuscript.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Drewes M. The phenomenology of right handed neutrinos. Int J Mod Phys E (2013) 22:1330019. doi:10.1142/s0218301313300191

CrossRef Full Text | Google Scholar

2. Atre A, Han T, Pascoli S, Zhang B. The search for heavy Majorana neutrinos. JHEP (2009) 05(030). doi:10.1088/1126-6708/2009/05/030

CrossRef Full Text | Google Scholar

3. King SF Neutrino mass models. Rept Prog Phys (2004) 67:107–57. doi:10.1088/0034-4885/67/2/r01

CrossRef Full Text | Google Scholar

4. Drewes M, Agostini M, Ky NA, Araki T, Archidiacono M, Bahr M, et al. A white paper on keV sterile neutrino dark matter. JCAP (2017) 2017:025. doi:10.1088/1475-7516/2017/01/025

CrossRef Full Text | Google Scholar

5. Davidson S, Nardi E, Nir Y. Leptogenesis. Phys Rept (2008) 466:105–77. doi:10.1016/j.physrep.2008.06.002

CrossRef Full Text | Google Scholar

6. Abazajian KN, Acero MA, Agarwalla SK, Albright CH, Antusch S, Balantekin AB, et al. Light sterile neutrinos: a white paper. High Energ Phys - Phenomenology (2012) 4. doi:10.48550/arXiv.1204.5379

CrossRef Full Text | Google Scholar

7. Berryman JM, de Gouvea A, Fox PJ, Kayser BJ, Kelly KJ, Raaf JL. Searches for decays of new particles in the DUNE multi-purpose near detector. JHEP (2020) 02:174. doi:10.1007/jhep02(2020)174

CrossRef Full Text | Google Scholar

8. Shrock RE. New tests for, and bounds on, neutrino masses and lepton mixing. Phys Lett B (1980) 96:159–64. doi:10.1016/0370-2693(80)90235-x

CrossRef Full Text | Google Scholar

9. Shrock RE. General theory of weak processes involving neutrinos. I. Leptonic pseudoscalar-meson decays, with associated tests for, and bounds on, neutrino masses and lepton mixing. Phys Rev D (1981) 24:1232–74. doi:10.1103/physrevd.24.1232

CrossRef Full Text | Google Scholar

10. Aguilar-Arevalo A, Aoki M, Blecher M, Britton D, vom Bruch D, Bryman D, et al. Search for heavy neutrinos in π → μν decay. Phys Lett B (2019) 798:134980. doi:10.1016/j.physletb.2019.134980

CrossRef Full Text | Google Scholar

11. Artamonov. Search for heavy neutrinos in K⁺ → μ⁺ν_H decays. Phys Rev D (2015) 91:052001. doi:10.1103/PhysRevD.91.052001

CrossRef Full Text | Google Scholar

12. Abi B, et al. Deep underground neutrino experiment (DUNE), far detector technical design report, volume I introduction to DUNE. JINST (2020) 15(08):T08008. doi:10.1088/1748-0221/15/08/T08008

CrossRef Full Text | Google Scholar

13. Abi B, et al. Deep underground neutrino experiment (DUNE), far detector technical design report. Vol. II DUNE Phys. (2020) 2. doi:10.48550/arXiv.2002.03005

CrossRef Full Text | Google Scholar

14. Giganti C, Lavignac S, Zito M. Neutrino oscillations: the rise of the PMNS paradigm. Prog Part Nucl Phys (2018) 98:1–54. doi:10.1016/j.ppnp.2017.10.001

CrossRef Full Text | Google Scholar

15. Gronau M, Leung CN, Rosner JL Extending limits on neutral heavy leptons. Phys Rev D (1984) 29:2539–58. doi:10.1103/physrevd.29.2539

CrossRef Full Text | Google Scholar

16. Bondarenko K, Boyarsky A, Gorbunov D, Ruchayskiy O. Phenomenology of GeV-scale heavy neutral leptons. JHEP (2018) 11:032. doi:10.1007/jhep11(2018)032

CrossRef Full Text | Google Scholar

17. Richman JD, Burchat PR Leptonic and semileptonic decays of charm and bottom hadrons. Rev Mod Phys (1995) 67:893–976. doi:10.1103/revmodphys.67.893

CrossRef Full Text | Google Scholar

18. Ballett P, Boschi T, Pascoli S. Heavy neutral leptons from low-scale seesaws at the DUNE near detector. JHEP (2020) 03:111. doi:10.1007/jhep03(2020)111

CrossRef Full Text | Google Scholar

19. Abada A, Bečirević D, Sumensari O, Weiland C, Zukanovich Funchal R Sterile neutrinos facing kaon physics experiments. Phys Rev D (2017) 95(7):075023. doi:10.1103/physrevd.95.075023

CrossRef Full Text | Google Scholar

20. Kayser B, Shrock RE Distinguishing between Dirac and Majorana neutrinos in neutral current reactions. Phys Lett B (1982) 112:137–42. doi:10.1016/0370-2693(82)90314-8

CrossRef Full Text | Google Scholar

21. Adam AA, Kordosky M, On behalf of the Dune Collaboration . Deep underground neutrino experiment (DUNE) near detector conceptual design report. Instruments (2021) 5(4):31. doi:10.3390/instruments5040031

CrossRef Full Text | Google Scholar

22. DUNE Collaboration. Home (2021). Available from: https://glaucus.crc.nd.edu/DUNEFluxes/.

23. Agostinelli S, Allison J, Amako K, Apostolakis J, Araujo H, Arce P, et al. GEANT4–a simulation toolkit. Nucl Instrum Meth A (2003) 506:250–303. doi:10.1016/s0168-9002(03)01368-8

CrossRef Full Text | Google Scholar

24. Allison J, Amako K, Apostolakis J, Arce P, Asai M, Aso T, et al. Recent developments in Geant4. Nucl Instrum Meth A (2016) 835:186–225. doi:10.1016/j.nima.2016.06.125

CrossRef Full Text | Google Scholar

25. Ferrari A, Paola R. Sala, Alberto fasso, and johannes ranft. FLUKA: a multi-particle transport code. Program version 2005 (2005) 10. doi:10.5170/CERN-2005-010

CrossRef Full Text | Google Scholar

26. Böhlen TT, Cerutti F, Chin MPW, Fassò A, Ferrari A, Ortega PG, et al. The FLUKA code: developments and challenges for high energy and medical applications. Nucl Data Sheets (2014) 120:211–4. doi:10.1016/j.nds.2014.07.049

CrossRef Full Text | Google Scholar

27. Sjöstrand T, Ask S, Christiansen JR, Corke R, Desai N, Ilten P, et al. An introduction to PYTHIA 8.2. Comput Phys Commun (2015) 191:159–77. doi:10.1016/j.cpc.2015.01.024

CrossRef Full Text | Google Scholar

28. Abi B, Acciarri R, Acero MA, Adamov G, Adams D, Adinolfi M, et al. Experiment simulation configurations approximating DUNE TDR. High Energ Phys - Exp (2021) 3. doi:10.48550/arXiv.2103.04797

CrossRef Full Text | Google Scholar

29. Bolton PD, Deppisch FF, Bhupal Dev PS Neutrinoless double beta decay versus other probes of heavy sterile neutrinos. JHEP (2020) 03:170. doi:10.1007/jhep03(2020)170

CrossRef Full Text | Google Scholar

30. Andreopoulos C, Bell A, Bhattacharya D, Cavanna F, Dobson J, Dytman S, et al. The GENIE neutrino Monte Carlo generator. Nucl Instrum Meth A (2010) 614:87–104. doi:10.1016/j.nima.2009.12.009

CrossRef Full Text | Google Scholar

31. Bischer I, Werner R. General neutrino interactions at the DUNE near detector. Phys Rev D (2019) 99(3):036006. doi:10.1103/physrevd.99.036006

CrossRef Full Text | Google Scholar

32. Bryman DA, Shrock R Improved constraints on sterile neutrinos in the MeV to GeV mass range. Phys Rev D (2019) 100(5):053006. doi:10.1103/physrevd.100.053006

CrossRef Full Text | Google Scholar

33. Argüelles CA, Foppiani N, Hostert M. Heavy neutral leptons below the kaon mass at hodoscopic neutrino detectors. Phys Rev D (2022) 105(9):095006. doi:10.1103/physrevd.105.095006

CrossRef Full Text | Google Scholar