Search for heavy neutrinos with the T2K near detector ND280

TL;DR

This study searches for heavy neutrinos in the 140-493 MeV/c^2 mass range using the T2K ND280 detector, setting new upper limits on their mixing with standard neutrinos by analyzing decay signatures and finding no evidence of their existence.

Contribution

It provides the first constraints on heavy neutrino mixing in this mass range using the T2K ND280 detector data, improving previous experimental limits.

Findings

No excess events observed in the signal region.

Set upper limits on heavy neutrino mixing elements, e.g., $U_e^2 < 10^{-9}$ at 90% C.L.

Constraints are competitive with prior experiments.

Abstract

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.

Tables

Table 1. Table 1: 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 subsection II.1 .

	Neutrino mode				Anti-neutrino mode
	CR-I		CR-II		CR-I		CR-II
Ch.	D	NEUT	D	NEUT	D	NEUT	D	NEUT
${\mu }^{\pm }{\pi }^{\mp }$	$15$	$11.4\pm 1.0$	$36$	$30.1\pm 1.6$	$2$	$2.6\pm 0.5$	$5$	$6.5\pm 0.8$
$e^{-}{\pi }^{+}$	$2$	$2.4\pm 0.5$	$14$	$13.3\pm 1.1$	$2$	$0.7\pm 0.2$	$7$	$5.4\pm 0.7$
$e^{+}{\pi }^{-}$	$5$	$3.1\pm 0.5$	$17$	$11.7\pm 1.0$	$1$	$0.5\pm 0.2$	$2$	$3.3\pm 0.6$
${\mu }^{+}{\mu }^{-}$	$0$	$0.9\pm 0.3$	$3$	$2.5\pm 0.5$	$0$	$0.2\pm 0.1$	$1$	$0.4\pm 0.2$
$e^{+}{e}^{-}$	$4$	$4.5\pm 0.6$	$58$	$53.5\pm 2.2$	$3$	$0.8\pm 0.3$	$13$	$14.3\pm 1.3$

Table 2. Table 2: Summary of the estimated numbers of background events in the different analysis channels in neutrino and anti-neutrino beam modes with the corresponding absolute uncertainties (statistical, flux-related, detector-related, model), using the data set presented in subsection II.1 .

Mode	Ch.	Expected	Uncertainties
		background	stat.	flux	det.	model	total
neutrino	${\mu }^{\pm }{\pi }^{\mp }$	1.543	0.366	0.154	0.165	0.285	0.516
	$e^{-}{\pi }^{+}$	0.376	0.213	0.038	0.104	0.097	0.259
	$e^{+}{\pi }^{-}$	0.328	0.186	0.033	0.117	0.115	0.250
	${\mu }^{+}{\mu }^{-}$	0.216	0.107	0.022	0.045	0.062	0.133
	$e^{+}{e}^{-}$	0.563	0.192	0.056	0.092	0.074	0.233
anti-neutrino	${\mu }^{\pm }{\pi }^{\mp }$	0.384	0.161	0.038	0.058	0.100	0.202
	$e^{-}{\pi }^{+}$	0.018	0.018	0.002	0.005	0.005	0.020
	$e^{+}{\pi }^{-}$	0.219	0.155	0.022	0.140	0.122	0.243
	${\mu }^{+}{\mu }^{-}$	0.038	0.038	0.004	0.007	0.011	0.040
	$e^{+}{e}^{-}$	0.015	0.015	0.002	0.001	0.004	0.016

Equations

$$\displaystyle\mathcal{L}_{\text{mass}}=-\frac{1}{2}\left(\begin{array}[]{cc}\bar{\nu}_{L}\ \bar{\nu}_{R}^{c}\end{array}\right)\left(\begin{array}[]{cc}0&m_{D}\\ m_{D}^{T}&m_{R}\end{array}\right)\left(\begin{array}[]{c}\nu_{L}^{c}\\ \nu_{R}\end{array}\right)+\text{h.c.},$$

$${\nu }_{\alpha }=\sum _{i=1}^3{V}_{\alpha i}{\nu }_i+\sum _{I=1}^n{\Theta }_{\alpha I}{N}_I(\alpha =e,\mu ,\tau ),$$

$${\mathcal{N}}_A=B_A+\sum _i{\epsilon }_{A,i}\times f_i(U_e^2,U_{\mu }^2,U_{\tau }^2)\times {\Phi }_i.$$

$$\mathcal{L}=\prod _A\text{Poisson}\left(n_A^{\text{obs}},{\mathcal{N}}_A\right).$$

$$p(U_e^2,U_{\mu }^2,U_{\tau }^2)=\int d\Phi d\epsilon dB\times \mathcal{L}\times {\pi }_{\Phi }{\pi }_{\epsilon }{\pi }_B{\pi }_{U^2}.$$

$$p_{\text{prof}}(U_e^2)$$

$$p_{\text{marg}}(U_e^2)$$

Full text

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 $\Delta m^{2}$) consists in the introduction of $n\geq 2$ new right-handed (sterile) neutrino fields $\nu_{R}$ with the following mass term Drewes (2013):

[TABLE]

where $m_{D}$ is the $3\times n$ Dirac mass matrix and $m_{R}$ is the $n\times n$ Majorana mass matrix. If the seesaw condition $m_{D}^{T}m_{D}\ll m_{R}^{2}$ holds (in terms of eigenvalues), diagonalisation of the mass matrix yields three light Majorana mass eigenstates $\nu_{i}$ ($i=1,2,3$), with masses $m_{\nu,i}$ of the order of the eigenvalues of $m_{D}m_{R}^{-1}m_{D}^{T}$ and $n$ heavy Majorana mass eigenstates $N_{I}$ ($I=1,\ldots,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:

[TABLE]

where $V$ corresponds to the usual PMNS matrix and $\Theta$ is the active-heavy mixing matrix. Heavy neutrinos can be produced in leptonic meson decays $M^{\pm}\to\ell_{\alpha}^{\pm}+N_{I}$ with a branching ratio proportional to $|\Theta_{\alpha I}|^{2}$, for $M_{N}<m_{\text{meson}}-m_{\ell_{\alpha}}$ and they can similarly decay via the same mixing element.

If at least two of the heavy neutrinos have a mass between $0.1$ and $100$ GeV/c2, they can generate baryogenesis via leptogenesis without any additional new physics Drewes et al. (2018). An example of such a model is the Neutrino Minimal Standard Model ($\nu$MSM) with $n=3$, in which $N_{1}$ has a mass of $1-100$ keV/c2 and is a warm dark matter candidate, while $N_{2,3}$ are degenerate with GeV-scale masses Asaka and Shaposhnikov (2005); Asaka et al. (2005).

In the following, we define $U_{\alpha}^{2}\equiv\sum|\Theta_{\alpha I}|^{2}$ summing over the heavy neutrinos that cannot be distinguished experimentally (such as $N_{2}$ and $N_{3}$ in the $\nu$MSM).

Limits on $U_{\alpha}^{2}$ for $M_{N}<493$ MeV/c2 can be obtained either by studying heavy neutrino production from kaon decays ($K^{\pm}\to\ell^{\pm}N$) or by searching for heavy neutrino decays, e.g. to one pion and one charged lepton ($N\to\ell_{\alpha}^{\pm}\pi^{\mp}$). The best constraints in this mass range were obtained by the BNL E949 Artamonov et al. (2015) and the CERN PS191 Bernardi et al. (1986, 1988) experiments with limits of the order of $10^{-9}-10^{-8}$ on $U_{e}^{2}$, $U_{\mu}^{2}$ and $U_{e}U_{\mu}$ for $M_{N}=200-450$ MeV/c2. Limits from other experiments are summarised in the review Drewes and Garbrecht (2017).

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 Asaka et al. (2013).

II Experimental setup

II.1 The T2K beamline

The Tokai-to-Kamioka (T2K) experiment Abe et al. (2011) 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\times 10^{20}$ protons-on-target (POT) in neutrino mode and $6.29\times 10^{20}$ POT in anti-neutrino mode, after data quality cuts.

II.2 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 Abe et al. (2011). The central tracker consists of three Time Projection Chambers (TPCs) Abgrall et al. (2011), two scintillator-based Fine-Grained Detectors (FGDs) Amaudruz et al. (2012) and one $\pi^{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 cross-section 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 Abe et al. (2015). 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$ m3.

III Analysis

III.1 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 Abe et al. (2013); Abgrall et al. (2016). We first consider 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^{\pm}\to\ell^{\pm}_{\alpha}N$, $\alpha=e,\mu$) by weighting event-by-event using the appropriate branching ratios Shrock (1981a, b); Gorbunov and Shaposhnikov (2007) and modified kinematics. The analysis assumes the heavy neutrino lifetime is long enough to reach ND280 ($\tau\gg\text{time of flight}\sim 1\,\mu$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 $\sim 3$ lower.

The heavy neutrino decays are then simulated at a random point along their trajectories inside ND280. All the possible modes $N\to\ell^{\pm}\pi^{\mp}$ and N\to\ell^{\pm}\ell^{\mp}\overset{\scalebox{0.3}{(}\raisebox{-1.7pt}[0.0pt][0.0pt]{-}\scalebox{0.3}{)}}{\nu} were simulated. Figure 3 shows the allowed production and decay modes as a function of the heavy neutrino mass. The neutral current decay modes N\to e^{+}e^{-}\overset{\scalebox{0.3}{(}\raisebox{-1.7pt}[0.0pt][0.0pt]{-}\scalebox{0.3}{)}}{\nu}_{\tau} and N\to\mu^{+}\mu^{-}\overset{\scalebox{0.3}{(}\raisebox{-1.7pt}[0.0pt][0.0pt]{-}\scalebox{0.3}{)}}{\nu}_{\tau} are directly sensitive to the mixing element $U_{\tau}^{2}$.

Effects related to heavy neutrino polarisation Levy (2018) and delayed arrival time (with respect to light neutrinos) are taken into account in the simulation.

III.2 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: $\mu^{\pm}\pi^{\mp}$, $e^{-}\pi^{+}$, $e^{+}\pi^{-}$, $e^{+}e^{-}$, $\mu^{+}\mu^{-}$.

In the analysis, we do not define any specific selection for the three-body decays $N\to e^{\pm}\mu^{\mp}\nu$ because these modes already contribute to the $e^{\pm}\pi^{\mp}$ selection channels. For the $\mu^{+}\mu^{-}$ 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_{\rm inv}$ of the two-track system: in the case of a heavy neutrino decay, it is expected that $m_{\rm inv}^{\rm true}\leq M_{N}$ ($m_{\rm inv}^{\rm true}=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/c2 allowing an upper cut on the reconstructed invariant mass $m_{\rm inv}^{\rm reco}<700$ MeV/c2 to be applied. The additional margin accounts for detector resolution effects.

• •

angle between the two tracks $\Delta\Phi$: the two charged tracks produced in the decay are boosted forward so that only events with $\Delta\Phi<90\degree$ can be selected without loss of signal efficiency.

• •

incoming heavy neutrino polar angle $\theta$: 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. $\theta$ 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\theta>0.992$ for $\mu^{\pm}\pi^{\mp}$ for the $\mu\pi$ channel and $\cos\theta>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^{\pm}\to\mu^{\pm}N$ or $K^{\pm}\to e^{\pm}N$). $\mu\pi$ efficiencies are slightly better as muon tracks are easier to reconstruct in the TPC.

III.3 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 subsection III.1, 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 Abe et al. (2013).

• •

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 Abe et al. (2018). The overall effects have been estimated to be approximately $5\%$.

III.4 Background estimation

The background remaining after the selection has first been estimated using the NEUT 5.3.2 Monte Carlo generator Hayato (2002), before being constrained using control regions in ND280 data.

One of the dominant background contributions is the neutrino-induced coherent pion production on argon nuclei in the TPC gas ($\nu_{\mu}+\text{Ar}\to\mu^{-}+\pi^{+}+\text{Ar}$). The NEUT prediction has been tuned to T2K and MINERvA data Abe et al. (2016); Higuera et al. (2014) 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 $\theta$ 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 1 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 2 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 $\mu\pi$ channel is higher than for other channels, as it is dominated by the irreducible coherent pion production.

III.5 Statistical analysis

Two approaches have been considered to constrain the mixing elements $U_{e}^{2}$, $U_{\mu}^{2}$ and $U_{\tau}^{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 $\mu^{\pm}\pi^{\mp}$ channel as defined in subsection III.2 can constrain:

• •

either $U_{\mu}^{2}$ by considering only the signal from $K^{\pm}\to\mu^{\pm}N$, $N\to\mu^{\pm}\pi^{\mp}$,

• •

or $U_{e}\times U_{\mu}$ by considering only the signal from $K^{\pm}\to e^{\pm}N$, $N\to\mu^{\pm}\pi^{\mp}$.

Three methods to obtain constraints in this approach have been applied:

1. ($\mathbf{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 Cousins and Highland (1992); 2. ($\mathbf{B}$)

the Feldman-Cousins method to define confidence intervals, taking into account the non-zero background Feldman and Cousins (1998); 3. ($\mathbf{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 Bernardi et al. (1986, 1988). 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_{e}^{2}$, $U_{\mu}^{2}$ and $U_{\tau}^{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_{e}^{2}=U_{\mu}^{2}=U_{\tau}^{2}=1$ and 100% selection efficiency, denoted $\Phi_{i}$ ;

• •

the selection efficiency of these decays in the current channel $\varepsilon_{A,i}$ ;

• •

the actual values of $U_{e}^{2}$, $U_{\mu}^{2}$ and $U_{\tau}^{2}$ via the factor $f_{i}=U_{\alpha}^{2}\sum U_{\beta_{j}}^{2}$ with $\alpha,\beta_{j}\in\{e,\mu,\tau\}$ where $\alpha$ is the flavour involved at the production of the heavy neutrino and $\beta_{j}$ are the flavours involved in its decay (only one for charge current modes, several for neutral current modes).

The expected number of events $\mathcal{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:

[TABLE]

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^{\text{obs}}_{A}$ in each channel $A$, with Poisson parameter $\mathcal{N}_{A}$:

[TABLE]

The uncertainties on the flux and efficiency are taken into account in the forms of multivariate Gaussian priors $\pi_{\Phi}$ and $\pi_{\varepsilon}$ respectively. The priors on the background $\pi_{B}$ are taken to be log-normal with means and standard deviations given by the expected background and its uncertainty in Table 2. The priors on the mixing elements $U_{\alpha}^{2}$ are assumed to be flat.

The marginalised posterior probability $p$ is then defined as the product of the likelihood $\mathcal{L}$ and the priors, integrating over all the nuisance parameters (flux, efficiency and background):

[TABLE]

A Markov Chain Monte Carlo method has been implemented using PyMC Patil et al. (2010) 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,

[TABLE]

where $U_{\mu,\text{max}}^{2}$ and $U_{\tau,\text{max}}^{2}$ are the values maximising the likelihood.

Limits in 2D/3D parameter space may be obtained as well. Limits on $U_{e}^{2}$ can be computed for $140<M_{N}<493$ MeV/c2, while limits on $U_{\mu}^{2}$ and $U_{\tau}^{2}$ can only be computed for $140<M_{N}<388$ MeV/c2 due to the kinematic constraints presented in Figure 3.

IV Results

Following the selection from subsection III.2, no events were observed in any of the different signal regions, which is consistent with the background-only hypothesis, allowing upper limits on $U_{e}^{2}$, $U_{\mu}^{2}$ and $U_{\tau}^{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 ($\mathbf{A}$, $\mathbf{B}$, $\mathbf{C}$), which give similar upper limits with method $\mathbf{A}$ giving slightly more conservative limits as expected.

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 Bernardi et al. (1986, 1988), E949 Artamonov et al. (2015) and CHARM Orloff et al. (2002), especially in the high-mass region (above $300$ MeV/c2). The kinks clearly visible on $U_{\mu}^{2}$ and $U_{\tau}^{2}$ limits come from the changes in the contributing production and decay modes as presented in Figure 3.

The limits are obtained after marginalisation over the two other mixing elements. For $U_{e}^{2}$, the limits after profiling (equation 12) are also presented, which effectively corresponds to setting $U_{\mu}^{2}=U_{\tau}^{2}=0$. Indeed, for $M_{N}>388$ MeV/c2, the correlations between $U_{e}^{2}$ and $U_{\mu}^{2}$ (as seen in Figure 7) would give limits on $U_{e}^{2}$ outside T2K’s reach. However, profiling leads to a loss in the sensitivity on $U_{e}^{2}$ with respect to the marginalisation as it forcefully suppresses the contributions of the decay modes involving $U_{\mu}^{2}$ or $U_{\tau}^{2}$.

It is worth mentioning that the limits depend on the choice of prior on $U_{\alpha}^{2}$. The limits on $U_{e}^{2}$ and $U_{\mu}^{2}$ are quite robust with respect to a change of prior as T2K data are directly sensitive to these mixing elements (e.g. using $\pi_{U^{2}}(U_{\alpha}^{2})=U_{\alpha}^{2}$ varies the limit by less than $30\%$), while the limit on $U_{\tau}^{2}$ is strongly affected (more than $50\%$).

It is also possible to define 2D contours, e.g. in the $U_{e}^{2}-U_{\mu}^{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/c2 corresponds to the kinematic cut-off for $K^{\pm}\to\mu^{\pm}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_{e}^{2}$, $U_{\mu}^{2}$ and $U_{\tau}^{2}$ are obtained using a combined Bayesian approach. Results apply to any model with heavy neutrinos with masses between $140$ MeV/c2 and $493$ MeV/c2 such as Rasmussen and Winter (2016), 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 $\nu$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/c2. 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/c2 and $U_{e,\mu}^{2}>10^{-11}-10^{-10}$) to be fully explored.

The results presented in this article are available in the corresponding data release dat (2018). 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.

Acknowledgments

We thank the J-PARC staff for superb accelerator performance. We thank the CERN NA61/SHINE Collaboration for providing valuable particle production data. We acknowledge the support of MEXT, Japan; NSERC (Grant No. SAPPJ-2014-00031), NRC and CFI, Canada; CEA and CNRS/IN2P3, France; DFG, Germany; INFN, Italy; National Science Centre (NCN) and Ministry of Science and Higher Education, Poland; RSF, RFBR, and MES, Russia; MINECO and ERDF funds, Spain; SNSF and SERI, Switzerland; STFC, UK; and DOE, USA. We also thank CERN for the UA1/NOMAD magnet, DESY for the HERA-B magnet mover system, NII for SINET4, the WestGrid and SciNet consortia in Compute Canada, and GridPP in the United Kingdom. In addition, participation of individual researchers and institutions has been further supported by funds from ERC (FP7), “la Caixa” Foundation (ID 100010434, fellowship code LCF/BQ/IN17/11620050), the European Union’s Horizon 2020 Research and Innovation programme under the Marie Sklodowska-Curie grant agreement no. 713673 and H2020 Grant No. RISE-GA644294-JENNIFER 2020; JSPS, Japan; Royal Society, UK; the Alfred P. Sloan Foundation and the DOE Early Career program, USA; RFBR research project 18-32-00072, Russia.