Search for a heavy neutral lepton that mixes predominantly with the tau neutrino

We report a search for a heavy neutral lepton (HNL) that mixes predominantly with ν τ . The search utilizes data collected with the Belle detector at the KEKB asymmetric energy e + e − collider. The data sample was collected at and just below the center-of-mass energies of the Υ (4 S ) and Υ (5 S ) resonances and has an integrated luminosity of 915 fb − 1 , corresponding to (836 ± 12) × 10 6 e + e − → τ + τ − events. We search for production of the HNL (denoted N ) in the decay τ − → π − N followed by its decay via N → µ + µ − ν τ . The search focuses on the parameter-space region in which the HNL is long-lived, so that the µ + µ − originate from a common vertex that is significantly displaced from the collision point of the KEKB beams. Consistent with the expected background yield, one event is observed in the data sample after application of all the event-selection criteria. We report limits on the mixing parameter of the HNL with the τ neutrino as a function of the HNL mass.

While the standard model (SM) has only left-handed, massless neutrinos, measurements of neutrino flavor oscillations [1] demonstrate that neutrinos do have mass.This suggests the possible existence of right-handed neutrino states, which carry no SM charges, with corresponding Yukawa and Majorana mass terms in the Lagrangian.A large ratio between the mass terms gives rise to a seesaw mechanism [2][3][4][5][6][7][8] that explains the small neutrino masses and predicts the existence of heavy mass eigenstates referred to as heavy neutral leptons (HNLs).In particular, discovery of an HNL in the GeV mass scale can also explain the origin of the baryon asymmetry of the Universe via leptogenesis [9][10][11] and would motivate searches for an additional, keV-scale HNL that could be a dark matter candidate [12][13][14][15][16].
Among the direct searches, both the BABAR [33] search and the current search utilize the large sample of e + e − → τ + τ − events available at e + e − B factories.The BABARsearch employs a missing-energy signature.The search reported here uses the displaced vertex (DV) formed by the decaying HNL, following the method of Ref. [37], which is used here for the first time.
We search for production of the HNL in the decay τ − → N π − .Being lighter than the τ lepton, the HNL undergoes the weak-neutral-current decay N → ν τ µ + µ − .Thus, the signal signature is a pion that is prompt, i.e., originating from a point near the e + e − interaction point (IP), and a µ + µ − pair originating from a DV, i.e., a point significantly displaced from the IP.The search is sensitive to a range of m N and |V N τ | 2 values for which the HNL is long-lived [38] yet decays within the Belle tracking system.
The Belle detector, described in detail in Refs.[39,40], is a large-solid-angle magnetic spectrometer that consists of a silicon vertex detector, a 50-layer central drift chamber (CDC) that spans the radial region 8.3 < r < 86.3 cm, an array of aerogel threshold Cherenkov counters (ACC), time-of-flight scintillation counters (TOF), and an electromagnetic calorimeter (ECL) composed of CsI(Tl) crystals located inside a superconducting solenoid coil that provides a 1.5 T magnetic field.Finally, the KLM system consists of the iron flux-return located outside of the coil, instrumented to detect K 0 L mesons and to identify muons.
Electron identification is performed with a normalized likelihood P e , determined from specific ionization (dE/dx) in the CDC, the number of photons in the ACC, the ECL cluster shape, the ratio of energy deposited in the ECL to the track momentum, and the track-ECL cluster position matching [43].Muon candidates are identified with a normalized likelihood P µ determined from dE/dx, ACC, TOF, and KLM information [44].
To optimize the event selection, model the background, and determine the signal efficiency, we use various Monte Carlo (MC) simulation samples.Background samples include the processes e + e − → τ + τ − , e + e − → q q where q is a u, d, s, or c quark, e + e − → B B, e + e − → e + e − , e + e − → µ + µ − , and e + e − → e + e − f f (where f is an electron, muon, or quark), generated with the event generators KKMC [45] with TAUOLA [46], PYTHIA [47], EvtGen [48] with PYTHIA, BHLUMI [49], KKMC, and AAFH [50], respectively.Most background-MC samples correspond to integrated luminosities that are 5 times that of the data.While the MC samples for e + e − → e + e − and e + e − → e + e − f f are smaller, the contribution of these events to the background is negligible.Signal-MC samples of e + e − → τ + τ − events are generated with KKMC.The generic decays of one τ lepton, referred to as the tag τ , are simulated with TAUOLA and PYTHIA according to the known branching fractions [1].The decays of the signal τ and of the HNL are simulated with PYTHIA.Signal MC samples are generated for m N values between 300 and 1600 MeV/c 2 in 25-MeV/c 2 -wide steps.Depending on m N , the HNL lifetime τ N is simulated such that cτ N is between 15 and 30 cm, chosen so that a large fraction of the simulated HNLs decay within the CDC.Photon production via initial-state radiation is simulated in all but the e + e − → B B events.In all MC samples, final-state radiation from charged particles is simulated with PHOTOS [51], and the detector-response simulation is performed with GEANT3 [52].
In the signal samples, the HNL decay is simulated with equal probability density in all regions of the µ + µ − ν τ phase space, ignoring spin correlations.To correct for this, signal events are also generated with Mad-Graph5 aMC@NLO [53] with the HeavyN model of Majorana neutrinos [54] and an effective vertex of the form ∂ µ πτ γ µ (1 − γ 5 )N + H.c for the τ − → π − N decay.For each value of m N , the PYTHIA-generated events are weighted according to the Dalitz-plot distribution of the MadGraph5 aMC@NLO events at the generator level.All other kinematic-variable distributions are consistent among the PYTHIA and MadGraph5 aMC@NLO events, and the impact on the efficiency of potential correlations between the angular distributions of the tag-τ and signal-τ decays are determined to be negligible.
The event-selection criteria are decided by comparing distributions of simulated signal and background events.Each event is divided into two hemispheres centered on the c.m. thrust axis [55,56], which is determined from the tracks and photons in the event.One hemisphere must contain exactly three tracks and the other exactly one track, taken as the signal-τ and tag-τ decays, respectively.The sum of the track charges must be zero.The prompt pion candidate and the tag-τ daughter-track candidate must have transverse and longitudinal impact parameters of |dr| < 0.5 and |dz| < 2.0 cm with respect to the IP.The pion candidate must satisfy P µ < 0.01 and P e < 0.01.Pairs of photons, each with energy above 0.1 GeV, are combined to form π 0 candidates if their invariant mass satisfies 0.115 < m(γγ) < 0.152 GeV/c 2 , corresponding to about ±3.5 standard deviations (σ) with respect to the m(γγ) resolution around the nominal π 0 mass [1].Events containing a π 0 candidate are discarded.The total lab-frame energy of photons with energies greater than 200 MeV that are not used to form π 0 candidates must be less than 1 GeV, to reject hadronic background while allowing for initial-state radiation, machine background, and detector noise in signal events.
HNL candidates are reconstructed from a µ + µ − pair, where the muons are selected with the requirement P µ > 0.9.The muons are fitted to a common vertex that constitutes the DV.Henceforth, the momenta of the muons are evaluated at the DV.The χ 2 probability of the DV fit is required to be greater than 10 −5 .The radial position of the DV must satisfy r DV > 15 cm with respect to the CDC symmetry axis, to suppress background from prompt tracks, decays of K 0 S mesons and hyperons, and particle interaction in dense detector material.To further suppress background from prompt tracks while retaining high signal efficiency, each muon must satisfy r hit − r DV > −2 cm, where r hit is the radial position of the smallest-radius CDC hit of the muon.The c.m.-frame angle between the muons must satisfy cos θ µ + µ − > 0.5 to suppress cosmic-ray background and ensure that the muons are consistent with originating from a boosted parent.To ensure consistency with the τ decay, we apply the invariant-mass requirement m πDV = (p π + p DV ) 2 /c < 1.776 GeV/c 2 , where p π and p DV are the 4-momenta of the pion and of the µ + µ − , respectively Events are rejected if they satisfy 420 < m DV (ππ) < 520 MeV/c 2 , where m DV (ππ) is the dimuon mass calculated with the pion mass hypothesis for the two muons.This rejects K 0 S → π + π − decays in which the pions decay to or are misidentified as muons.
Despite the neutrino in the final state, the constraints of the signal decay enable reconstruction of the full kinematics of the signal-τ decay chain with a twofold ambiguity [37], which arises from a quadratic equation in the magnitude p N of the HNL momentum, Here we have defined (with all quantities in the laboratory frame) where m τ , m π , and m 2 µµ are the masses of the τ , π ± [1], and µµ system, E = E π + E µµ , q X = p X cos θ N X , p X and E X are the measured momentum and energy for X = π or X = µµ, and θ N X is the angle between pX and pN .We take pN to be the vector between the IP and the DV, neglecting the small τ flight distance.The derivation of these expressions is explained in the Appendix [57].The discriminant S of Eq. ( 1) tends to be larger for background than for signal.We require it to satisfy S < 0.4 GeV 2 , retaining more than 98% of signal events.Given a solution for p N , the squared HNL mass is The two solutions for m N are referred to as m + and m − , depending on the sign in front of the square root of the quadratic-equation solution.For events with S < 0 we set S = 0, in which case m + = m − .Events that pass these selections are divided into two signal regions (SRs).The region SRH, defined by the requirement m DV (ππ) > 520 MeV/c 2 , targets heavy HNLs.Light HNLs are targeted by the region SRL, which is defined by m DV (ππ) < 420 MeV/c 2 .Furthermore, events in the SRL are required to satisfy either m + < 900 or m − < 600 MeV/c 2 .
In the background-MC sample, we find two background events in the SRH, both from e + e − → τ + τ − .One event contains the decay τ − → ν τ π − K + K − , and the DV is formed from the muons originating from the decays π − → µ − νµ and K + → µ + ν µ .The other event contains τ − → ν τ π − K 0 S , with both pion daughters of the K 0 S decaying to the muons that constitute the DV.In the SRL the MC background yield is two e + e − → τ + τ − events and two e + e − → q q events.In all four events, the DV is formed from a muon produced in a pion decay and either another such muon or a pion from a K 0 L or K 0 S decay.Given that the integrated luminosity of the data is 1/5 that of the MC samples, the MC prediction for the background yield is 0.40 ± 0.28 in the SRH and 0.80±0.40 in the SRL, where the uncertainties arise from the finite MC sample size.
In addition to the use of MC samples, we study the background from the event yields in control regions (CRs) labeled CRH and CRL, in correspondence to the signal regions SRH and SRL.Each CR is defined by the same selection criteria as the corresponding SR, except that only one of the tracks forming the DV is identified as a muon with P µ > 0.9, while the other must satisfy P e < 0.1 and P µ < 0.1 to suppress leptons and enhance the pion contribution.In the background-MC sample, the CRH has 337 e + e − → τ + τ − events.The DVs in these events are formed from pions produced in K 0 S → π + π − decays (in 76% of the events), K 0 S → π + π − followed by pion decay to a muon (13%), prompt π ± or K ± mesons (2%) or muons produced in their decays (4%), K 0 L decays (2%) and hadronic interactions in material (2%).The CRH also contains 26 e + e − → q q events, where the DV is formed from K 0 S decay (38%), Λ → pπ − (27%), π ± or K ± decays (12%), particles produced in hadronic interaction in detector material (12%), K 0 L decay (8%), and prompt π ± or K ± (4%).The CRL has 101 e + e − → τ + τ − , in which the DVs are formed from K 0 L decay (57%), e + e − from photon conversion (26%), prompt π ± or K ± (15%) and K 0 L decay (2%).The CRL also has 81 e + e − → q q events, with DVs formed from Λ → pπ − (59%), K 0 L → πµν (40%), and K 0 S decays (1%).Several validation regions (VRs), titled VRHππ, VRLππ, VRK S , VRHss, and VRLss, are used for data-MC consistency studies and systematic uncertainty estimation.The VRHππ and VRLππ are defined by the same selection criteria as the SRH and SRL, respectively, except that both DV tracks must satisfy the nonlepton requirement P e < 0.1, P µ < 0.1.Data-MC comparison in these regions helps validate the background estimated for DVs produced from pairs of hadrons.The VRK S , defined identically to VRHππ but with 480 < m DV (ππ) < 515 GeV/c 2 , is used to check the overall level of background from K S decays.The VRHss and VRLss are defined by the same criteria as SRH and SRL, respectively, except that the electric charges of the two muons have the same sign, opposite the charge of the prompt pion.These regions are used to validate the level of potential background from coincidental crossing of tracks.The event types that populate the VRs in the MC are listed in the Appendix [57].
For each region, the event yields observed in the data, the corresponding MC prediction, their ratio, and the statistical consistency between them are shown in Table I.Also shown are the postfit yields in the SRs and CRs for the case of no signal.The data contain one event in the SRH, with m + = m − = 1.473GeV/c 2 , and no events in the SRL.To avoid potential experimenter bias, the data event yields in the SRs were unveiled only after finalizing all analysis procedures and systematic uncertainty estimations.In the CRs and VRs, the number of data events exceed the MC expectation by between 2% and 43%, and the naive data-MC statistical consistency N σ obs,bgd = (N obs − N bgd )/ N obs + σ 2 bgd ranges between 0.8 and 4.7.While decays of τ leptons are well simulated, MC-data differences may arise from the simulation of q q events, which is not as well tested for low-multiplicity events.The largest of these differences is used to estimate an uncertainty on the background model, described later.
The S values of the data and background-MC events that pass all the selection criteria except S < 0.4 are shown in Fig. 1(a) together with the signal-MC distribution.The m ± values after all selections except the m ± requirement are shown in Fig. 1(b).The SRL distribution for signal MC events with m N = 600 MeV/c 2 is also shown for comparison.Signal events cluster around either m + ≈ m N or m − ≈ m N , with events for which S was set to 0 having m − = m + .
From the observed event yields N obs in the SRs and CRs we compute 95% confidence-level (CL) upper limits on |V N τ | 2 as a function of m N using the CL s prescription implemented in pyhf [58][59][60] with the likelihood function The expected event yield in region sig , where N R bgd is the expected background yield shown in Table I, and N R sig is the expected signal yield, calculated as Here the index i indicates data samples with different values of E c.m. ; L i is the integrated luminosity for sample i; σ i is the cross section for e + e − → τ + τ − , taken to be 0.919 ± 0.003 nb for E c.m. = E c.m. (Υ (4S)) and scaled by E 2 c.m. (Υ (4S))/E 2 c.m. for other samples [61]; B(τ → πN ) and B(N → µµν) are the branching fractions of the specified decays [38] listed in the Appendix [57]; and ϵ R is the total efficiency for signal events to satisfy all the selection criteria of region R, determined from the signal-MC samples.Since muons have a finite probability of failing the muon-identification criteria, signal events may populate the CRs.The ratio N CR sig /N SR sig between the signal yields in each CR and its corresponding SR ranges between 0.1 and 2.1, depending on the HNL mass.Nonetheless, the  The m− vs. m+ values for the data and for e + e − → τ + τ − and e + e − → q q background-MC events (which has 5 times the data luminosity) in the SRH (black symbols) and SRL (red symbols).The region enclosed by the dashed lines is vetoed in the SRL.The distribution of signal-MC events with mN = 600 MeV/c 2 in the SRL is also shown in colored contours.
CRs are dominated by background events for all values of |V N τ | 2 for which the likelihood is not negligible.
The efficiency ϵ R depends on m N and the HNL lifetime, which is taken from Ref. [38] and listed in the Appendix [57] for each value of m N and |V N τ | 2 .While the signal-MC events are generated with a specific lifetime value τ 0 N , we determine ϵ R for any lifetime τ N , as follows.First, we use each signal-MC sample to calculate the efficiency as a function of the radial and longitudinal position (r DV , z DV ) of the HNL decay, creating an efficiency map ϵ R (m N ; r DV , z DV ).The map bin size is 5 cm in z DV and 10 cm in r DV .Larger bins are used for the SRL (SRH) efficiency of heavy (light) HNLs, for which the sample size is limited.For each event in the signal-MC sample we randomly draw a set of decay times t i from the exponential distribution exp(−t i /τ N )/τ N .For each t i we calculate the event's HNL decay position and determine its efficiency from the efficiency map.The individual event efficiencies are used to determine the total signal efficiency given the MadGraph5 aMC@NLO weights.The resulting signal efficiencies in the SRs are given in the Appendix.[57].The calculation of ϵ R , N R sig , and the upper limit is performed on a grid in |V N τ | 2vs-m N space.The grid uses a 25 MeV/c 2 step in m N , corresponding to the generated signal-MC samples.In |V N τ | 2 the grid is logarithmic, with a multiplicative step size of 10 0.05 .
All systematic uncertainties are handled with the nuisance parameters shown in Eq. ( 4).The largest relative systematic uncertainty, 34%, is assigned to the background yield expectations N R bgd .This value is chosen since it brings the event yields in the data and MC samples, N VRHππ obs and N VRHππ exp , to within 1σ consistency in the VRHππ, which has the poorest data-MC consistency, as shown in Table I.Following Ref. [17], we assign a 5% uncertainty on the HNL branching fraction and decay modeling, arising mainly from the QCD corrections to the HNL hadronic decay width [38,62].The systematic uncertainty on the integrated luminosity is 1.4% [63], and that on the e + e − → τ + τ − cross section is 0.3% [61].The uncertainty on the reconstruction of the two prompt tracks is 0.7% [64].The relative statistical uncertainties on ϵ R associated with the finite number of MC events, the finite number of generated decay times t i , and the MadGraph5 aMC@NLO weights are also used as systematic uncertainties.The uncertainty due to muon identification is 2% per muon [64].The uncertainty on the efficiency for the remaining event-reconstruction steps, namely, online event selection and trigger, and tracking and vertexing of the HNL daughter tracks, is estimated to be 3.7%, chosen so that it brings the data and MC yields to within 1σ consistency in the VRK S .We estimate an uncertainty of 1.3% on the determination of the signal efficiency from the maximal difference with respect to the true efficiency at the generated lifetime.The postfit values of the nuisance parameters are all well within 1σ of their expected values, and are shown in the Appendix [57].The largest pull is that of the backgroundprediction parameter µ B , 0.69 ± 0.43.
The resulting excluded region in |V N τ | 2 -vs-m N space are shown in Fig. 2, separately for a Dirac and a Majorana HNL.For every point on the grid, the HNL lifetime in the Majorana case is half that of the Dirac HNL.
In conclusion, we report a search for a heavy neutral lepton in the decay chain The search method, used here for the first time, utilizes the displaced vertex originating from the long-lived HNL decay and the ability to reconstruct the HNL-candidate mass to suppress the background to the single-event level.It also allows for direct measurement of the HNL mass if a signal is observed.In the signal regions targeting heavy and light HNLs we observe 1 and 0 events, respectively, FIG.2: The expected (dashed) and observed (solid) 95% CL limits on |VNτ | 2 vs. mN for a Dirac or Majorana HNL.The green and yellow bands show the 1σ and 2σ bands for the expected limits for the Dirac case.The blue and pink bands show the same for the Majorana case.Also shown are the limits from direct searches at DELPHI [26], corrected for the unavailability of the charged-current decay for mN < mτ [65], ArgoNeuT [32] (90% CL), and the upper limit from BABAR [33].The 90% CL limits arising from reinterpretations [34,35] of other searches by CHARM [29,36] and WA66 [30] are shown as well.
in agreement with the background expectation.We set limits on the mixing coefficient |V N τ | 2 of the HNL with the τ neutrino for HNL masses in the range 300 < m N < 1600 MeV/c 2 .
Appendix: Efficiency. Figure 3 shows the total efficiency for signal selection as a function of the mixing parameter |V N τ | 2 and the mass m N .Examples of the efficiency maps are shown for m N = 300, 1600 MeV/c 2 in Fig. 4. The efficiency for the various off-line selection requirements (cut flow) is shown in Table II for selected signal samples.Origin of DV tracks in MC events in the VRππ regions.Table III gives details of the origin of the DV tracks in e + e − → τ + τ − and e + e − → q q MC events.
Estimation of Nuisance Parameter pulls.Figure 5 shows the pulls for each of the nuisance parameters when the expected signal yield is zero.
Calculation of the HNL mass values m ± .Consider the decay chain τ → N x, N → ν τ y, where the kinematics of the x and y systems are measured.We start with 4momentum conservation in the τ decay, p τ = p N + p x .Squaring this relates the HNL energy E N , momentum |⃗ p N |, and mass m N to those of the x and to the τ mass, where Here θ N x is the angle between the momentum directions of the x and the HNL, the latter given by the position   showing the selection efficiency at different stages of the off-line selection for three signal samples, with masses of 300, 1000, and 1600 MeV/c 2 , generated with cτ values of 15, 22.5, and 30 cm, respectively.The first row gives the efficiency for the online selection, which includes the requirement r DV > 5 cm, loose muon identification for the muon candidates, and that the tag-τ track is in the opposite hemisphere to the three signal-τ tracks.In the leftmost column, the final selections for SRH and SRL are labeled.

Region
VRHππ VRLππ Event sample τ τ q q τ τ q q Total no.  of the DV relative to the IP, where the small flight distance of the τ can be neglected.Similarly, 4-momentum conservation in the HNL decay yields where and θ N y is the angle between the momenta of the y the HNL.Equations ( 6) and ( 8) give E N in terms of |⃗ p N |: where are known.The HNL mass is then and is also given by inserting Eq. ( 10) in ( 6), where are known.Comparing Eqs. ( 12) and ( 13) gives the quadratic equation with the solution where S = (2AB − D) 2 − 4(B 2 − 1)(A 2 − C).Inserting Eq. ( 16) into Eq.( 13) gives the two solutions for m N .HNL lifetime and branching fractions.Table IV lists the values used for the HNL lifetime, its production branching fraction in τ decay, and the branching fraction for its decay into µµν.
Here the index R runs through the regions SRH, SRL, CRH, and CRL; P N R obs |N R exp is the Poisson probability for an observation of N R obs events in region R given the expectation N R exp ; the index C runs through the nuisance parameters; and G C (p C |p 0 C , σ C ) is the Gaussian distribution for nuisance parameter p C given the expectation p 0 C and its uncertainty σ C .

FIG. 1 :
FIG. 1: (a)The S values of the data and MC events after applying all SR requirements except S < 0.4.The signal-MC distribution is arbitrarily normalized, and the background-MC distributions are normalized to the data luminosity.(b) The m− vs. m+ values for the data and for e + e − → τ + τ − and e + e − → q q background-MC events (which has 5 times the data luminosity) in the SRH (black symbols) and SRL (red symbols).The region enclosed by the dashed lines is vetoed in the SRL.The distribution of signal-MC events with mN = 600 MeV/c 2 in the SRL is also shown in colored contours.

FIG. 3 :
FIG. 3: The total signal efficiencies in (top) SRH and (bottom) SRL as a function of |VNτ | 2 and mN .

FIG. 4 :
FIG. 4: Final reconstruction efficiency for signal events as a function of the radial and longitudinal positions of the DV for (top) mN = 1600 MeV/c 2 events in the SRH and (bottom) mN = 300 MeV/c 2 events in the SRL.

FIG. 5 :
FIG.5:The pull values of the nuisance parameters for the fit to data with µ fixed to 0. The nuisance parameters correspond to the systematic uncertainties arising from (shown from left to right) the background prediction in SRH and SRL; the product of the integrated luminosity L, e + e − → τ + τ − cross section σ, and the branching fractions Bτ and BN for the signal τ and HNL decays; and the background sample sizes in the CRH, CRL, SRH, SRL.

TABLE I :
The number N obs of events observed in the data, the expected number N bgd of background events based on the MC simulation, the ratio N obs /N bgd , and the naive statistical consistency N σ obs,bgd = (N obs − N bgd )/ N obs + σ 2 bgd for each of the signal, control, and validation regions.The last column shows the postfit event yields for the SRs and CRs, obtained by maximizing the likelihood function when the signal yield is fixed to 0.

TABLE II :
Cut flow table

TABLE III :
Origin of the tracks forming the DV in the VRHππ and VRLππ validation regions in the MC sample, whose integrated luminosity is 5 times that of the data.The different categories are mutually exclusive.In the "hadronic material interaction" and K/π decay in flight categories either one or both tracks may be of this origin.