New constraints on tau-coupled Heavy Neutral Leptons with masses $m_N = 280-970$ MeV

A search for Heavy Neutral Leptons has been performed with the ArgoNeuT detector exposed to the NuMI neutrino beam at Fermilab. We search for the decay signature $N \to \nu \mu^+ \mu^-$, considering decays occurring both inside ArgoNeuT and in the upstream cavern. In the data, corresponding to an exposure to $1.25 \times 10^{20}$ POT, zero passing events are observed consistent with the expected background. This measurement leads to a new constraint at 90\% confidence level on the mixing angle $\left\vert U_{\tau N}\right\rvert^2$ of tau-coupled Dirac Heavy Neutral Leptons with masses $m_N =$ 280 - 970 MeV, assuming $\left\vert U_{eN}\right\rvert^2 = \left\vert U_{\mu N}\right\rvert^2 = 0$.

Introduction.-The discovery that neutrinos oscillate and therefore have mass has inspired numerous experimental efforts to understand this phenomenon. The Standard Model (SM) does not predict the existence of neutrino masses, requiring additional fields and/or interactions to generate them. One such model requires the existence of two or more Heavy Neutral Leptons (HNLs): SM gauge singlet fermions that mix with the light neutrinos. This mixing can induce the observed small neutrino masses via one of many different seesaw mechanisms [1][2][3][4][5][6][7]. In addition, HNLs can provide solutions to other mysteries of nature such as the baryon asymmetry of the universe [8] (via leptogenesis) or dark matter [9]. In this paper, we present a search for HNLs with masses O(100) MeV using the ArgoNeuT detector.
We consider the simplest phenomenological scenario including a HNL, N -that it has a mass m N and mixes with the light neutrinos via one or more non-zero new angles |U eN | 2 , U µN 2 , and |U τ N | 2 in an extended 4 × 4 leptonic mixing matrix. If m N is in the ∼MeV-GeV range, HNLs can be produced as a result of high-energy protonfixed-target collisions, travel to a downstream detector and decay producing detectable charged particles. In the ArgoNeuT detector, we search for the decay signature N → νµ + µ − .
ArgoNeuT was a 0.24 ton Liquid Argon Time Projection Chamber (LArTPC) neutrino detector located in the NuMI beam [10] at Fermilab that collected data in [2009][2010]. The instrumented volume of the TPC was 40×47×90 cm 3 (vertical, drift, beam direction) with two readout planes, each consisting of 240 wires spaced by 4 mm and oriented at ±60 • to the horizontal. A detailed description of the design and operation of the ArgoNeuT detector can be found in Ref. [11]. The ArgoNeuT detector was located 100 m underground in the MINOS near detector hall, 1033 m downstream of the NuMI target and immediately upstream of the MINOS near detector (MINOS-ND). ArgoNeuT was able to use the MINOS-ND as a muon spectrometer. A detailed description of the MINOS-ND can be found in Ref. [12]. The analysis reported in this paper is performed using 1.25 × 10 20 protons-on-target (POT) collected in reverse horn current (anti-neutrino) mode, during which both ArgoNeuT and the MINOS-ND were operational [11].
Generation and simulation.-A HNL with mass m N and mixing angle |U αN | 2 with the light neutrinos, ν α (α = e, µ, τ ), can be produced by any kinematically accessible process that would normally result in an outgoing ν α . For the decay N → νµ + µ − we require m N > 2m µ and, for simplicity, we assume that only one angle |U αN | 2 is non-zero at a time. A variety of experiments [13][14][15][16][17][18][19][20][21] have set powerful constraints on the angles |U eN | 2 and U µN 2 in the region of interest for ArgoNeuT. We therefore focus on the case where |U τ N | 2 is the only non-zero mixing angle. In this scenario, the HNLs are predominantly produced in the decays of τ ± leptons originating arXiv:2106.13684v2 [hep-ex] 17 Aug 2021 from decays of D ± (s) mesons. The lifetime of N , as well as the branching ratio Br(N → νµ + µ − ), can be calculated as a function of the mixing |U τ N | 2 considering all kinematically accessible final states [22]. We assume that N is a Dirac fermion throughout this analysis.
In the NuMI beam approximately 87% of the incident 120 GeV protons interact in the target, with the majority of the remaining 13% interacting 715 m downstream in the hadron absorber [10]. We consider HNL production occurring in both the target and the absorber, the latter giving access to shorter N lifetimes as a result of being significantly closer to the detector. We simulate the particle propagation using GEANT4 [23] and the τ ± production using PYTHIA8 [24]. Approximately 10% of the beam protons reach the absorber with energy T p ≈ 120 GeV. For 120 GeV protons interacting in either the target or the absorber, an average of 2.1 × 10 −7 (3.0 × 10 −7 ) τ + (τ − ) are produced per proton [25]. To generate a flux of N , we simulate the decays τ ± → N X, where X consists of SM particles. We simulate the kinematics by assuming m X = m π ± and that the branching ratio of this new decay is Br(τ ± → N X ± ) = 0.9 |U τ N | 2 K(m N ) [22,26]. Since the D ± (s) and τ ± lifetimes are small, the kinematics of N produced in the target and absorber are qualitatively the same. However, the geometric acceptance of ArgoNeuT is significantly larger for the absorber-produced N due to the proximity to the detector.
The HNL decay products are then simulated in the ArgoNeuT detector using the LArSoft software framework [27], which simulates the particle propagation using GEANT4 [23] then performs detector response simulation and reconstruction [11,28]. A stand-alone version of the MINOS simulation and reconstruction is then used to simulate the tracks exiting ArgoNeuT and entering the MINOS-ND.
Signature.-The HNL decay N → νµ + µ − is seen in ArgoNeuT as a pair of minimally ionising particles (MIPs) that can be matched to a pair of oppositely charged particles in the MINOS-ND. These muons are energetic and highly forward-going: with average energy E µ ± ∼ 7 GeV; average angle with respect to the beam direction θ beam ∼ 1.5 • ; and an average opening angle θ opening ∼ 3 • . Given the ArgoNeuT angular resolution of approximately 3 • [29], this results in the muon pair frequently overlapping and being reconstructed as a single track for part or all of their length. Two in-ArgoNeuT decay signatures are therefore considered, each of which is illustrated in Fig. 1 (top, middle). In the first, the muons are reconstructed as two distinct MIP tracks originating from a common vertex, each of which can be matched to tracks in the MINOS-ND. This signature will be referred to as a two-track event. In the second, the muons overlap and are reconstructed as a single track with double-MIP dE/dx for part or all of their length. Then, in the MINOS-ND, the pair of oppositely-charged muons separate due to the presence of a magnetic field. This signature will be referred to as a double-MIP event.
In addition to decays occurring inside the ArgoNeuT detector, we also consider decays occurring in the cavern upstream of ArgoNeuT along the NuMI beam-line where the resulting muons then pass through the detector. This scenario is illustrated in Fig. 1 (bottom). During the ArgoNeuT physics run, the MINERvA detector [30] was under construction in the upstream cavern. We therefore only consider decays that occur in the 63 cm between the end of the MINERvA detector and the start of the ArgoNeuT TPC. In this scenario, only the double-MIP signature is considered. This is because the two-track signature is more difficult to distinguish from neutrinoinduced background muons due to the absence of vertex information, whereas the double-MIP signature is unique to potential HNL decays.
An example simulated HNL decay with a double-MIP signature is shown in Fig. 2. The pair of muon tracks fully overlap in ArgoNeuT (top, middle), then split once reaching the MINOS-ND (bottom) due to the magnetic field. The strongest identifier of whether a pair of overlapping muons are present is provided by the dE/dx of the track. The region of interest is the start of each track, prior to the muon pair possibly splitting. Fig. 3 shows the average reconstructed dE/dx over the first 5 cm of tracks resulting from simulated HNL decays. Two distinct peaks are visible. The first is at dE/dx ∼ 2 MeV/cm, approx- imately the dE/dx of a single minimally-ionizing muon. For these events the opening angle of the muons is sufficiently large to properly reconstruct them as two separate tracks. The second peak is at dE/dx ∼ 4.5 MeV/cm, approximately double the single MIP dE/dx, indicating two overlapping muons. A threshold is applied between the two peaks, illustrated by the dashed line, separating the double-MIP-like and MIP-like populations.
Selection.-A series of pre-selection cuts are first applied to remove poorly reconstructed events along with obvious non-HNL interactions. The highly forward-going muons from HNL decays can be challenging to reconstruct correctly in LArTPC detectors. This is because the ionisation tracks are near parallel to the readout planes and hence the drifted ionisation charge arrives on the wires at approximately the same time. These events may have large regions missed during the reconstruction and cannot be reliably identified. Therefore, we first remove events with fewer than 80% of reconstructed energy depositions associated with reconstructed tracks. Next, events with more than three total tracks with length L ≥ 5 cm, or more than two tracks originating from a vertex are removed. This removes any events that are obvious non-HNL interactions due to having additional reconstructed particles present. A harsher cut requiring only two tracks to be present is not applied because a common failure mode of the reconstruction is the presence of split tracks in the region where the overlapping muons begin to separate. Tracks shorter than L = 5 cm are not considered to avoid removing events that have short δ-rays originating from the muons. Events passing the pre-selection are then assessed against the two-track and double-MIP selection criteria sequentially.
In the two-track scenario only tracks starting within a fiducial volume in ArgoNeuT are considered, defined as: 1 ≤ x ≤ 46 cm (drift), −19 ≤ y ≤ 19 cm (vertical) and z ≥ 3 cm (beam direction), to remove backgrounds originating from the cavern. Events with two tracks that either originate from or can be projected back to a common vertex within the fiducial volume are selected. The tracks are required to be forward-going with respect to the beam direction, have length L ≥ 5 cm, exit Ar-goNeuT towards the MINOS-ND, have a mean dE/dx over their full length consistent with being a single MIP (dE/dx < 3.1 MeV/cm) and have an opening angle between them of θ opening ≤ 10 • .
In the double-MIP scenario, HNL decays occurring both inside ArgoNeuT and in the upstream cavern are considered. Any events containing tracks with an angle with respect to the beam direction θ beam > 15 are removed, as these likely originate from background interactions. The average dE/dx is then calculated over the first 10 hits (∼ 5 cm) of each track, where any individual anomalously large hits (dE/dx > 10 MeV/cm) are discarded. A cut is applied at dE/dx > 3.1 MeV/cm, illustrated by the dashed line in Fig. 3, to identify events with a possible pair of overlapping muons.
Once candidate events are identified with either the two-track or double-MIP signature in ArgoNeuT, MINOS-ND matching is performed. Each track is projected to the start of the MINOS-ND and the radial and angular off-sets between the projected tracks and each reconstructed MINOS-ND track are compared. In the two-track case, ArgoNeuT-MINOS-ND matching tolerances of r dif f ≤ 12.0 cm and θ dif f ≤ 0.17 rad are used [28]. In the double-MIP case, since a single track is being matched to two tracks in the MINOS-ND, the matching tolerances are loosened to 2.5 times the two-track case. The matched tracks are required to be forwardgoing with respect to the beam direction, start within 20 cm of the up-stream face of the detector and within the calorimeter region, and be at least 1 m long. This helps to remove any tracks that are unlikely to have originated from ArgoNeuT.
Finally, several selection cuts are applied in the MINOS-ND. These cuts are the same for both the twotrack and double-MIP scenarios. We require that the tracks have an average dE/dx consistent with being a muon (4 ≤ dE/dx ≤ 18 MeV/cm), are reconstructed with opposite charges, and have start times, t 0 , consistent with having originated from the same interaction or decay: |∆t 0 | ≤ 20 ns. Pairs of track with larger ∆t 0 could not have originated from a single HNL decay and instead are likely neutrino-induced background muons.
The selection efficiency as a function of the HNL energy, E N , is shown in Fig. 4 for simulated m N = 450 MeV HNL decays occurring inside the ArgoNeuT detector and at two positions in the upstream cavern. The efficiency inside the detector, defined as the fraction of events that are selected with either the two-track or double-MIP signatures, is around 60-65% and relatively flat above E N ∼ 10 GeV. However, it drops significantly at lower energies predominantly due to one or both of the muons being too low energy to reach the MINOS-ND. The cavern efficiencies are defined as the fraction of decays resulting in muons intersecting with the ArgoNeuT detector that are selected with the double-MIP signature. The further away from ArgoNeuT the decay occurs, the less likely the muon pair is to remain overlapping. This probability decreases further at lower energies where the muons are less forward-going.
Backgrounds and systematic uncertainties.-The primary backgrounds in this search originate from misreconstructed neutrino interactions occurring within the ArgoNeuT cryostat, and from neutrino-induced throughgoing muons arising from interactions upstream of the detector. Simulation of these backgrounds is performed with the GENIE [31] neutrino event generator using NuMI beam fluxes provided by the MINERvA collaboration [32], along with a data-driven model of neutrinoinduced through-going muons [29,33,34]. In the two track scenario, the dominant form of observed background events are charged current ν µ interactions where either the interaction vertex or one or more tracks have been poorly reconstructed leading to these events not being removed. We expect to see 0.1 ± 0.1 events of this type in the data. In the double-MIP scenario, the domi-  Table I. They are dominated by the uncertainty on the HNL flux. There is a 20% uncertainty on the D ± (s) production [35][36][37]. Then, the uncertainty on the branching ratios D ± (s) → τ ± + ν τ [38] leads to an additional 5.7% uncertainty on the τ ± flux. Combining these in quadrature leads to a 20.8% uncertainty on the resulting HNL flux. Next, we consider the impact of uncertainties in the reconstruction by repeating the analysis with each parameter varied individually according to its assigned uncertainty. We apply uncertainties of 3% on the tuning of the calorimetry [28], 3% on the track angular reconstruction [29] and 6% on the energy reconstruction of stopping particles in the MINOS-ND [39]. Finally, we assign a 1% uncertainty on the charge reconstruction due to the modelling of the magnetic field [39]. Combining the impact of the performed variations in quadrature leads to a 0.5% systematic uncertainty due to reconstruction effects. In addition to the reconstruction uncertainties, a 3.3% systematic uncertainty is assigned to the selection efficiency to account for the potential impact of neutrino-induced through-going muons present in 3.3% of triggers [11,40]. A through-going muon registered in coincidence with a HNL event would lead to it being discarded in the pre-selection. Finally, there is a 2.2% uncertainty in the size of the ArgoNeuT  [42] and DELPHI [43] are also shown in purple and blue, respectively.
instrumented volume originating from uncertainty in the electron drift velocity [29] and a 1% uncertainty in the number of collected POT [34].
Results.-The selection has been applied to the full ArgoNeuT 1.25 × 10 20 POT anti-neutrino mode data-set. In total zero events pass, consistent with the expected background rate of 0.4 ± 0.2 events. Figure 5 shows our exclusion of parameter space at 90% confidence level with 1.25×10 20 POT at ArgoNeuT, assuming production from τ ± decays. The limit is evaluated using a Bayesian approach with a uniform prior [41]. The ±1σ uncertainty on the expected constraint includes both the uncertainty on the background expectation and the 21.2% systematic uncertainty on the signal production, combined conservatively. The existing limits from CHARM [42] and DEL-PHI [43] are also shown in purple and blue, respectively. Our result leads to a significant increase in the exclusion region on the mixing angle |U τ N | 2 of tau-coupled Dirac HNLs with masses m N = 280 -970 MeV, assuming |U eN | 2 = U µN 2 = 0. Other scenarios are considered in the Supplemental Material [44].
Conclusions.-We have presented the first search for HNLs decaying with the signature N → νµ + µ − in a LArTPC detector. Applying a novel technique to iden-tify pairs of overlapping highly forward-going muons, we have searched for tau-coupled Dirac HNLs produced in the NuMI beam and decaying in the ArgoNeuT detector or in the upstream cavern. In the data, corresponding to an exposure to 1.25 × 10 20 POT, zero passing events are observed consistent with the expected background. The results of this search lead to a significant increase in the exclusion region on the mixing angle |U τ N | 2 of taucoupled Dirac HNLs with masses m N = 280 -970 MeV, assuming |U eN | 2 = U µN 2 = 0. The analysis techniques we developed could be applied in future HNL searches performed in larger mass LArTPC experiments.

ALTERNATIVE MODEL SCENARIOS
Here we consider two alternative scenarios where the only non-zero mixing angle between the Heavy Neutral Lepton, N , and the light neutrinos, ν α (α = e, µ, τ ), is either |U eN | 2 or |U µN | 2 . In each case we include two production channels, K ± → ± α N and D ± (s) → ± α N and, as before, the production rate of N is proportional to the mixing angle squared |U αN | 2 . Details about these production mechanisms can be found in Ref. [1]. The K ± channel is simulated using a GEANT4 [2] based simulation of the NuMI beam (G4NuMI), allowing for decay-in-flight of the focused K ± . In this case, N production is considered from interactions in the target only. The promptly decaying D ± (s) channel is simulated using PYTHIA8 [3] as in the main text, considering production in both the target and the absorber. The signal considered in these scenarios is the same as in the main text, N → νµ + µ − . However, the accepted HNL flux from these production mechanisms is different due to the different energy distribution of the parent particles that the HNLs originate from. This is accounted for using the energy-dependent efficiencies for the HNL selection discussed in the main text with the different energy distributions.
The systematic uncertainties on the HNL flux in the electron-and muon-coupled scenarios are different to the tau-coupled scenario discussed in the main text. In the case of HNLs originating from K ± decays, the uncertainty on the K ± production is +1.9 −5.8 % [4]. In the case of HNLs originating from D ± (s) decays, the uncertainty on the D ± (s) production remains the same, but the uncertainty on the resulting HNL flux is no longer impacted by the uncertainty on the D ± (s) → τ ± + ν τ branching ratio. The remaining systematic uncertainties are the same as discussed in the main text. Figure 1 shows our exclusion of parameter space at 90% confidence level in the scenarios where |U eN | 2 (left) and |U µN | 2 (right) are the only non-zero mixing angles. The uncertainties on the expected constraints at ±1σ assuming zero observed events are shown by the dashed black lines. In each case the strongest existing constraints are shown. In the electron-coupled scenario, these are from NA62 (red) [5], PS191 (green) [6,7], T2K (orange) [8], and CHARM (purple) [9]. In the muon-coupled scenario, these are from E949 (purple) [10], NA62 (red) [11], and NuTeV (blue) [12]. In both cases, less sensitive constraints are not shown to improve clarity. We find that the extracted limits from ArgoNeuT are considerably weaker than other existing limits. This is due to two factors: greater exposure and access to other final-states with larger branching ratios, such as N → e − π + in the electron-coupled case.