Search for heavy neutrinos with the T2K near detector ND280

This paper reports on the search for heavy neutrinos with masses in the range $140<M_N<493$ MeV/c$^2$ using the off-axis near detector ND280 of the T2K experiment. These particles can be produced from kaon decays in the standard neutrino beam and then subsequently decay in ND280. The decay modes under consideration are $N \to \ell^{\pm}_{\alpha} \pi^{\mp}$ and $N \to \ell^+_{\alpha} \ell^-_{\beta} \nu (\bar\nu)$ ($\alpha,\beta=e,\mu$). A search for such events has been made using the Time Projection Chambers of ND280, where the background has been reduced to less than two events in the current dataset in all channels. No excess has been observed in the signal region. A combined Bayesian statistical approach has been applied to extract upper limits on the mixing elements of heavy neutrinos to electron-, muon- and tau- flavoured currents ($U_e^2$, $U_{\mu}^2$, $U_{\tau}^2$) as a function of the heavy neutrino mass, e.g. $U_e^2<10^{-9}$ at $90\%$ C.L. for a mass of $390$ MeV/c$^2$. These constraints are competitive with previous experiments.

This paper reports on the search for heavy neutrinos with masses in the range 140 < MN < 493 MeV/c 2 using the off-axis near detector ND280 of the T2K experiment. These particles can be produced from kaon decays in the standard neutrino beam and then subsequently decay in ND280. The decay modes under consideration are N → ℓ ± α π ∓ and N → ℓ + α ℓ − β ( − ) ν (α, β = e, µ). A search for such events has been made using the Time Projection Chambers of ND280, where the background has been reduced to less than two events in the current dataset in all channels. No excess has been observed in the signal region. A combined Bayesian statistical approach has been applied to extract upper limits on the mixing elements of heavy neutrinos to electron-, muon-and tau-flavoured currents (U 2 e , U 2 µ , U 2 τ ) as a function of the heavy neutrino mass, e.g. U 2 e < 10 −9 at 90% C.L. for a mass of 390 MeV/c 2 . These constraints are competitive with previous experiments.

I. INTRODUCTION
Neutrino oscillations provide strong evidence that neutrinos are massive particles. Although in the minimal Standard Model they are massless, the most natural extension to allow non-zero masses compatible with oscillation experiments results (two different ∆m 2 ) consists in the introduction of n ≥ 2 new right-handed (sterile) neutrino fields ν R with the following mass term [1]: (1) where m D is the 3 × n Dirac mass matrix and m R is the n × n Majorana mass matrix. If the seesaw condition m T D m D ≪ m 2 R holds (in terms of eigenvalues), diagonalisation of the mass matrix yields three light Majorana mass eigenstates ν i (i = 1, 2, 3), with masses m ν,i of the order of the eigenvalues of m D m −1 R m T D and n heavy Majorana mass eigenstates N I (I = 1, . . . , n) (heavy neutrinos, also called heavy neutral leptons in the literature), with masses M N,I of the order of the eigenvalues of m R . The flavour eigenstates can be expressed in terms of the mass eigenstates as: Θ αI N I (α = e, µ, τ ), (2) where V corresponds to the usual PMNS matrix and Θ is the active-heavy mixing matrix. Heavy neutrinos can be produced in leptonic meson decays M ± → ℓ ± α + N I with a branching ratio proportional to |Θ αI | 2 , for M N < m meson − m ℓα and they can similarly decay via the same mixing element. * now at CERN † also at J-PARC, Tokai, Japan ‡ affiliated member at Kavli IPMU (WPI), the University of Tokyo, Japan § also at National Research Nuclear University "MEPhI" and Moscow Institute of Physics and Technology, Moscow, Russia ¶ deceased * * also at JINR, Dubna, Russia † † also at BMCC/CUNY, Science Department, New York, New York, U.S.A.
If at least two of the heavy neutrinos have a mass between 0.1 and 100 GeV/c 2 , they can generate baryogenesis via leptogenesis without any additional new physics [2]. An example of such a model is the Neutrino Minimal Standard Model (νMSM) with n = 3, in which N 1 has a mass of 1 − 100 keV/c 2 and is a warm dark matter candidate, while N 2,3 are degenerate with GeV-scale masses [3,4].
In the following, we define U 2 α ≡ |Θ αI | 2 summing over the heavy neutrinos that cannot be distinguished experimentally (such as N 2 and N 3 in the νMSM).
Limits on U 2 α for M N < 493 MeV/c 2 can be obtained either by studying heavy neutrino production from kaon decays (K ± → ℓ ± N ) or by searching for heavy neutrino decays, e.g. to one pion and one charged lepton (N → ℓ ± α π ∓ ). The best constraints in this mass range were obtained by the BNL E949 [5] and the CERN PS191 [6,7] experiments with limits of the order of 10 −9 − 10 −8 on U 2 e , U 2 µ and U e U µ for M N = 200 − 450 MeV/c 2 . Limits from other experiments are summarised in the review [8].
This paper presents the search for potential heavy neutrinos produced in the T2K decay volume and decaying in the T2K Near Detector, ND280, as was originally suggested in [9].

A. The T2K beamline
The Tokai-to-Kamioka (T2K) experiment [10] is a long-baseline neutrino experiment located in Japan with the primary goal of measuring muon (anti-)neutrino oscillations using Super-Kamiokande as its far detector. The T2K neutrino beam is produced at the Japan Proton Accelerator Research Complex (J-PARC) by colliding 30 GeV protons on a graphite target. The pions and kaons produced are focused and selected by charge with magnetic horns and subsequently decay in flight to neutrinos. Depending on the polarity of the current in the horns, the experiment can be run either in neutrino or anti-neutrino mode.
In this analysis, the production of heavy neutrinos from kaon decays in data taken from November 2010 to May 2017 are considered. This corresponds to a total exposure of 12.34 × 10 20 protons-on-target (POT) in neu-trino mode and 6.29 × 10 20 POT in anti-neutrino mode, after data quality cuts.

B. The off-axis near detector ND280
The off-axis near detector ND280 is located 280 metres from the proton target. It is composed of several sub-detectors with a 0.2 T magnet [10]. The central tracker consists of three Time Projection Chambers (TPCs) [11], two scintillator-based Fine-Grained Detectors (FGDs) [12] and one π 0 detector (P0D). It is surrounded by an Electromagnetic Calorimeter (ECal) and a Side Muon Range Detector (SMRD). A schematic view of ND280 is shown in Figure 1.
The main goal of ND280 is to detect neutrino interactions in order to constrain both neutrino flux and crosssection parameters. The TPCs are filled with a gas mixture based on argon gas and provide excellent track and momentum reconstruction with a typical resolution of 8% for 1 GeV/c tracks [13]. This can be combined with energy loss (dE/dx) measurements in order to perform particle identification (PID) of charged tracks crossing the TPCs. The analysis focuses on heavy neutrino decays occurring in the ND280 TPC gas volumes, which corresponds to a total volume of interest of 6.3 m 3 .

A. Simulation
The simulation of heavy neutrino production and decay is performed using the T2K neutrino flux prediction, which is constrained by the NA61/SHINE experiment results and by in-situ measurements [14,15]. We first con-sider the flux of standard light neutrinos coming from kaon decays in the beamline and crossing the ND280 TPCs. This flux is transformed into a flux of heavy neutrinos (K ± → ℓ ± α N , α = e, µ) by weighting eventby-event using the appropriate branching ratios [16][17][18] and modified kinematics. The analysis assumes the heavy neutrino lifetime is long enough to reach ND280 (τ ≫ time of flight ∼ 1 µs), which is consistent with current limits on the mixing elements. Figure 2 presents the results of the simulation for different heavy neutrino masses and for both production modes in neutrino mode. The flux has the same shape for anti-neutrino mode, although it is a factor of ∼ 3 lower. TPCs from K ± → µ ± N and K ± → e ± N for several values of MN , with the T2K beam in neutrino-mode and for U 2 e = U 2 µ = 1. The black dotted curve corresponds to the limiting case of a massless neutrino (N = ν); the one for K → eν is not drawn as it is a few orders of magnitude lower due to helicity suppression.
The heavy neutrino decays are then simulated at a random point along their trajectories inside ND280. All the possible modes N → ℓ ± π ∓ and N → ℓ ± ℓ ∓ ( − ) ν were simulated. Figure 3 shows the allowed production and decay modes as a function of the heavy neutrino mass. The neutral current decay modes ν τ are directly sensitive to the mixing element U 2 τ . Effects related to heavy neutrino polarisation [19] and delayed arrival time (with respect to light neutrinos) are taken into account in the simulation.

B. Selection
The selection was developed to isolate the signal events listed in Figure 3 from the background expected from standard neutrino interactions with matter. In order to significantly improve the signal to background ratio, which is inversely proportional to the density of the medium, only events occurring in the TPC gas volume are considered for this analysis.
Events are pre-selected by identifying two tracks of opposite charge originating from a vertex in a TPC. There should be no other tracks in the TPC itself or in the detector located directly upstream (e.g. P0D for the first TPC or the first FGD for the second TPC). Particle identification for each individual track is performed using energy loss in the TPC. Five channels are then identified: In the analysis, we do not define any specific selection for the three-body decays N → e ± µ ∓ ν because these modes already contribute to the e ± π ∓ selection channels. For the µ + µ − channel, electromagnetic calorimeter information is also used to clearly identify the two muons.
Several kinematic cuts are then applied to further reject the background: • invariant mass m inv of the two-track system: in the case of a heavy neutrino decay, it is expected that m true inv ≤ M N (m true inv = M N for the two-body decays). The heavy neutrino is produced in kaon decays so that it is necessarily lighter than M K = 493 MeV/c 2 allowing an upper cut on the reconstructed invariant mass m reco inv < 700 MeV/c 2 to be applied. The additional margin accounts for detector resolution effects.
• angle between the two tracks ∆Φ: the two charged tracks produced in the decay are boosted forward so that only events with ∆Φ < 90°can be selected without loss of signal efficiency.
• incoming heavy neutrino polar angle θ: the heavy neutrino's direction is collinear to the beam, while the products of an active neutrino interaction are expected to be distributed with a larger angle because of potential missed tracks or nuclear effects. θ is reconstructed using the two charged tracks but can still be used with a good approximation for the three-body decays. The cut is cos θ > 0.992 for µ ± π ∓ for the µπ channel and cos θ > 0.99 for the others.
Applying these criteria to the signal simulated in the ND280 TPC gas volumes, the efficiencies of the signal selection for the different modes were obtained, as shown in Figure 4. For a given mass, they are quite independent of the production mode (K ± → µ ± N or K ± → e ± N ). µπ efficiencies are slightly better as muon tracks are easier to reconstruct in the TPC.

C. Signal systematic uncertainties
Two sources of systematic uncertainties on the heavy neutrino signal are considered: • flux: uncertainties on the kaon flux used as input to the simulation, as presented in section III A, are directly transposed into uncertainties on the flux of heavy neutrinos reaching ND280. The total normalisation uncertainty has been estimated to be 15%, using external data such as those from the NA61/SHINE experiment [14].
• signal selection efficiency: detector systematic uncertainties are defined to cope with any discrepancies between data and simulation of the detector effects. The dominant uncertainties are related to TPC reconstruction and particle identification performances and have been computed as in previous ND280 analyses [20]. The overall effects have been estimated to be approximately 5%.

D. Background estimation
The background remaining after the selection has first been estimated using the NEUT 5.3.2 Monte Carlo gen- I. Comparison of number of events in data (D) and corresponding NEUT prediction (with statistical uncertainties) in the control regions used to determine the model uncertainties in the different channels, using the data set presented in section II A.
One of the dominant background contributions is the neutrino-induced coherent pion production on argon nuclei in the TPC gas (ν µ + Ar → µ − + π + + Ar). The NEUT prediction has been tuned to T2K and MINERvA data [22,23] with a 30% normalisation uncertainty.
Additional background sources include other types of neutrino interactions in the gas and interactions outside the gas. An example of the latter is the conversion of a photon, emitted by a neutrino interaction in a FGD, to an electron-positron pair.
Data and simulations are compared with two sets of control regions in order to estimate the model uncertainty on the background. First, a selection of events similar to the signal events, but where the kinematic cut on the polar angle θ is inverted (CR-I), contains mostly resonant pion production and quasi-elastic processes on argon. Similarly, control regions are identified by considering events starting in the borders of the TPC (meaning the box containing the gas) as the volume of interest rather than in the gas itself (CR-II). These control regions are dominated by photon conversions and other mis-reconstructed processes. Table I presents the comparison of T2K data and NEUT predictions in the aforementioned control regions. We have not found any significant discrepancies between data and Monte Carlo predictions in any of them. Conservatively, we have assigned to each background source a model uncertainty equal to the statistical uncertainty of the data in the corresponding control region.
For a given channel, the number of expected background events is the nominal value from NEUT and the total uncertainty is the sum of the contributions from the Monte Carlo statistical uncertainty, the flux and detector systematic uncertainties and the model uncertainties described above.
Table II summarises the background in the different analysis channels. The dominant contribution to its uncertainty comes from the limited statistics of the samples. The background in the µπ channel is higher than for other channels, as it is dominated by the irreducible coherent pion production.

E. Statistical analysis
Two approaches have been considered to constrain the mixing elements U 2 e , U 2 µ and U 2 τ . In the first approach, each heavy neutrino production/decay mode is considered independently and the corresponding analysis channel is used to put limits on the associated mixing elements. For instance, the µ ± π ∓ channel as defined in section III B can constrain: • either U 2 µ by considering only the signal from • or U e × U µ by considering only the signal from Three methods to obtain constraints in this approach have been applied: (A) assuming that the background is zero, set conservative upper limits, independently of background modelling and estimation, on the mixing elements using the Highland-Cousin method [24]; (B) the Feldman-Cousins method to define confidence intervals, taking into account the non-zero background [25]; (C) a Bayesian method to define credible intervals, taking into account the non-zero background.
This "single-channel" approach has the advantage of being straightforward and is similar to that of the PS191 collaboration [6,7]. However, it does not allow the different modes and channels to be combined, so that the constraints are valid only under strong assumptions of the hierarchy of U 2 e , U 2 µ and U 2 τ . A "combined" approach was then defined, in which all the heavy neutrino production and decay modes (presented in Figure 3) and the ten different analysis channels (five for each beam mode) are considered simultaneously. For a given analysis channel A, the contribution of a mode i is characterised by: • the expected number of decays in the detector assuming U 2 e = U 2 µ = U 2 τ = 1 and 100% selection efficiency, denoted Φ i ; • the selection efficiency of these decays in the current channel ε A,i ; • the actual values of U 2 e , U 2 µ and U 2 τ via the factor βj with α, β j ∈ {e, µ, τ } where α is the flavour involved at the production of the heavy neutrino and β j are the flavours involved in its decay (only one for charge current modes, several for neutral current modes).
The expected number of events N A in a channel A depends on the background in this channel B A and the sum of the contributions from the different production and decay modes: Only a Bayesian method has been considered in this combined approach. The likelihood is built using a Poisson function for the observed number of events n obs A in each channel A, with Poisson parameter N A : The uncertainties on the flux and efficiency are taken into account in the forms of multivariate Gaussian priors π Φ and π ε respectively. The priors on the background π B are taken to be log-normal with means and standard deviations given by the expected background and its uncertainty in Table II. The priors on the mixing elements U 2 α are assumed to be flat. The marginalised posterior probability p is then defined as the product of the likelihood L and the priors, integrating over all the nuisance parameters (flux, efficiency and background): A Markov Chain Monte Carlo method has been implemented using PyMC [26] to perform this integration. The output can then be used to define 90% domains, either by profiling or by marginalising over the two other mixing elements. For instance, where U 2 µ,max and U 2 τ,max are the values maximising the likelihood.
Limits in 2D/3D parameter space may be obtained as well. Limits on U 2 e can be computed for 140 < M N < 493 MeV/c 2 , while limits on U 2 µ and U 2 τ can only be computed for 140 < M N < 388 MeV/c 2 due to the kinematic constraints presented in Figure 3.

IV. RESULTS
Following the selection from section III B, no events were observed in any of the different signal regions, which is consistent with the background-only hypothesis, allowing upper limits on U 2 e , U 2 µ and U 2 τ to be placed. An example of results from the single-channel approach is presented in Figure 5. It shows the comparison of the three methods (A, B, C), which give similar upper limits with method A giving slightly more conservative limits as expected. e as a function of heavy neutrino mass using the single-channel approach, considering only the contribution from K ± → e ± N, N → e ± π ∓ , with the three methods A, B and C. The limits are compared to the ones of PS191 experiment [6,7].
The results of the combined approach are shown in Figure 6. They provide an improvement by a factor of 2-3 with respect to the single-channel approach, thanks to the increased statistical power of the combination.
The limits are competitive with those of previous experiments such as PS191 [6,7], E949 [5] and CHARM [27], especially in the high-mass region (above 300 MeV/c 2 ). The kinks clearly visible on U 2 µ and U 2 τ limits come from the changes in the contributing production and decay modes as presented in Figure 3. µ (middle), U 2 τ (bottom) as a function of heavy neutrino mass, obtained with the combined approach. The blue dashed lines corresponds to the results of the single-channel approach (method C). The blue solid lines are obtained after marginalisation over the two other mixing elements. In the top plot, the additional blue dotted line corresponds to the case where profiling is used (U 2 µ = U 2 τ = 0), as explained in the main text. The limits are compared to the ones of other experiments: PS191 [6,7], E949 [5], CHARM [27].
The limits are obtained after marginalisation over the two other mixing elements. For U 2 e , the limits after profiling (equation 7) are also presented, which effectively corresponds to setting U 2 µ = U 2 τ = 0. Indeed, for M N > 388 MeV/c 2 , the correlations between U 2 e and U 2 µ (as seen in Figure 7) would give limits on U 2 e outside T2K's reach. However, profiling leads to a loss in the sensitivity on U 2 e with respect to the marginalisation as it forcefully suppresses the contributions of the decay modes involving U 2 µ or U 2 τ . It is worth mentioning that the limits depend on the choice of prior on U 2 α . The limits on U 2 e and U 2 µ are quite robust with respect to a change of prior as T2K data are directly sensitive to these mixing elements (e.g. using π U 2 (U 2 α ) = U 2 α varies the limit by less than 30%), while the limit on U 2 τ is strongly affected (more than 50%). It is also possible to define 2D contours, e.g. in the U 2 e −U 2 µ plane, allowing the correlations between the mixing elements to be visualised. Figure 7 presents a set of such contours for different heavy neutrino masses. The change of behaviour at M N = 388 MeV/c 2 corresponds to the kinematic cut-off for K ± → µ ± N processes as seen in Figure 3.

V. CONCLUSION
A selection of events with two tracks with opposite charges originating from the ND280 TPC gas volumes allows heavy neutrino decays to be efficiently isolated from expected background coming from standard neutrino interactions with matter. No events are observed in the defined signal regions, which is consistent with the background-only hypothesis.
Limits on the mixing elements U 2 e , U 2 µ and U 2 τ are obtained using a combined Bayesian approach. Results apply to any model with heavy neutrinos with masses between 140 MeV/c 2 and 493 MeV/c 2 such as [28], and can, in particular, be interpreted as constraints on the sum of N 2 and N 3 coupling squared as explained in the introduction, for the νMSM.
As the analysis is still statistically limited, results are expected to further improve by a factor of 2-3 with T2K data up to 2026. Additional data will also allow the background treatment to be improved by using more populated control regions.
By considering heavy neutrino production from pion decays, it would also be possible to extend the phase space down to a few MeV/c 2 . When combined with a better understanding of the expected background, this may permit the low-mass heavy neutrino phase space (10 < M N < 493 MeV/c 2 and U 2 e,µ > 10 −11 − 10 −10 ) to be fully explored.
The results presented in this article are available in the corresponding data release [29]. It contains the signal flux and selection efficiencies for all modes and masses, the detailed background predictions, the limits presented in Figures 5-6-7 and the raw output of the MCMC. One can use the latter to re-compute the limits with different priors on the mixing elements or with different ways to present the results.