Constraints from the CHARM experiment on Heavy Neutral Leptons with tau mixing

We re-analyze the results of the searches for Heavy Neutral Leptons (HNLs) by the CHARM experiment. We study HNL decay channel $N\to e^{+}e^{-}\nu/\mu^{+}\mu^{-}\nu$ and show that, in addition to the constraints on the HNL's mixings with $\nu_e$ or $\nu_{\mu}$, the same data also implies limits on the HNLs that mix only with $\nu_{\tau}$ and have masses in the range $380 \text{ MeV}<m_{N}\lesssim 1.6\text{ GeV}$ - the region in the parameter space that was considered in the literature as a target for HNLs searches.


I. INTRODUCTION AND SUMMARY
The Standard Model (SM) of particle physics, despite its unprecedented success, has proven to be incomplete, encouraging searches for new particles.Heavy Neutral Leptons, or HNLs, is a well-motivated extension of the SM that is capable of simultaneously explaining several longstanding problems: neutrino oscillations [1], dark matter [2][3][4], and baryon asymmetry of the Universe [5][6][7][8].From the phenomenological point of view, such particles participate in weak interactions and behave as heavy neutrinos with interaction strength suppressed by the mixing angles U α , as compared to ordinary neutrinos ν α (see, e.g.[9]).HNLs with the masses in the MeV-GeV range are searched for at the accelerators (see, e.g.[10,11]), and may be constrained from cosmological observations as well [12][13][14].In the minimal models where two or three HNLs explain neutrino flavor oscillation data via a seesaw mechanism, each of the HNLs should have comparable mixings with all active neutrino flavors (see, e.g.[15][16][17] and references therein for discussion).Nevertheless, to obtain the accelerator bounds, it is convenient to consider HNLs that mix with only one flavor of active neutrinos, namely ν α (below, we denote such HNLs as N α ).Limits on U 2 α derived for such simplified "pure mixing model" are conservative, as the presence of additional mixing angles U 2 β =α = 0 for the same particle would only increase the number of expected events in a given experiment at the lower bound of sensitivity.The current bounds on U 2 e and U 2 τ for HNLs with pure mixing are shown in Fig. 2 as reported in [11].In the GeV mass range, the constraints on the mixing angle U 2 τ are orders of magnitude weaker as compared to the constraints on U 2 e/µ (constraints for the µ mixing are similar to the ones for the e mixing).Namely, for the e/µ mixing, the large values of the couplings for HNLs boiarska@nbi.ku.dk boyarsky@lorentz.leidenuniv.nlmikulenko@lorentz.leidenuniv.nlovchynnikov@lorentz.leidenuniv.nlwith masses m K m N m D 2 GeV are excluded by the CHARM experiment [18,19], while for the τ mixing CHARM constraints on U τ are reported in the literature only for masses m N < 290 MeV.The reason is the following: the original analysis [18,19] is based on negative results for searches for decays of feebly interacting particles into one of the possible dilepton pair -µe, µµ, µe.For HNLs, they consider only decays mediated through the charged current (CC) interaction (see Fig. 1, diagram (a)) that give rise to leptonic decays If only CC interactions are taken into account, the search is suitable to constrain the mixing of HNLs with ν e and ν µ .To search for CC mediated decays via the τ mixing (which necessarily include a τ lepton), the HNL mass should be m N > m τ m D in this model.Such HNLs are mainly produced in decays of heavy B mesons, the number of which at CHARM is insufficient to provide enough events for the couplings that are not excluded (see Fig. 2).Therefore, HNLs that mix only with ν τ cannot be constrained by CHARM data using only the decays via CC.In order to constrain the τ mixing angles of the light HNLs m N < m τ , one should include the interactions via the neutral current (NC) into the analysis, see Fig. 1 (diagram (b)).In this case, the dileptonic decays are and do not require the creation of a τ lepton for the pure τ mixing.
The works [16,20,21] have re-analyzed the CHARM constraints on HNLs by including also the neutral current processes.However, their analysis was insufficient to put the bounds on the pure τ mixing in GeV mass range.Namely, the work [20] (the results of which are used in [11]) has limited the study of the mass range by m N < 290 MeV, while [16,21] considered the decays of HNLs via neutral currents but did not include the production of HNLs from τ lepton (the diagrams (c) and (d) in Fig. 4).As a result, these works did not report any CHARM limits on the pure τ mixing.[19,20], NA62 [22], T2K [23], Belle [24,25], DELPHI [26], NOMAD [27], ArgoNeuT [28], see also [11] for a review and the references therein.For the pure τ mixing, we do not show the constraints imposed by the T2K experiment (unlike it is done in [11]), since they are reported for non-zero couplings U e/µ which dominate the production a .Constraints from the CHARM experiment are taken from the literature, while our re-analysis for them is shown in Fig. 6.The light gray domain corresponds to couplings that are either excluded by BBN [12,14] or too small to provide active neutrino masses.For the pure τ mixing, we also show sensitivities of the next generation Intensity frontier experiments (see text for details).In cyan, we show HNL parameter space that may be probed by neutrino observatories: the solid line shows the sensitivity of IceCube to the "double bang" signature from [29], while the dashed line corresponds to the sensitivity of KM3NeT to decays of HNLs produced in the atmosphere, see text and Appendix B for details.
a The T2K experiment [23] is based on a search for decays of HNLs produced from kaons K → lαN .Such decays can occur only through the e/µ mixings due to the small mass of kaon m K = 493 MeV.
To conclude, in the GeV mass range of HNLs, there is a gap in the parameter space probed by past experiments, possibly only due to the lack of the analysis (see Fig. 2).To close this gap, the searches for HNLs that mix mainly with τ neutrinos were considered among scientific goals for several experiments: displaced decays at FASER [10,30], Belle II [25], SND@LHC [31], DarkQuest [32], and NA62 in the dump mode [10]; prompt decays at LHCb [33,34]; and double bang signature at IceCube, SuperKamiokande, DUNE and Hyper-Kamiokande [29,35].
The planned neutrino observatory KM3NeT [36] working as an atmospheric beam dump may have sensitivity to such HNLs as well.Namely, HNLs may be produced in numerous collisions of cosmic protons with atmospheric particles, then reach the detector volume located deeply underwater in the Mediterranean Sea, and further decay into a dimuon pair inside.Such combination of decay products may be in principle distinguished from the SM events due to neutrino scatterings and penetrating atmospheric muons.We discuss this signature in more detail and estimate the sensitivity of KM3NeT to HNL produced in the atmosphere in Appendix B, and make the conclusions in Sec.IV.The resulting sensitivity is shown in Fig. 2. Yet, the main result of this paper is to remove the abovementioned limitations of the analysis for HNLs interacting only with τ neutrinos existing in the literature, and demonstrate that the results of the CHARM experiment do imply the constraints on the τ -mixing of HNLs with the masses in the range 290 MeV < m N < 1.6 GeV that are two orders of magnitude stronger than previously reported in the literature.Our results are shown in Fig. 6.
The CHARM bounds re-analysis presented in this paper may by similarly applied for the re-analysis of bounds coming from the NOMAD experiment [27].However, due to the smaller intensity of the proton beam at NOMAD and simultaneously similar geometric acceptance of the decay volume, the bounds imposed by NOMAD are subdominant, and we therefore do not make the re-analysis in this work.

II. CHARM EXPERIMENT
The CHARM experiment [18,19] was a proton beam dump operating at the 400 GeV CERN SPS.The total number of exposed protons was split into 1.7 • 10 18 protons on a solid copper target and 0.7•10 18 on a laminated copper target with the 1/3 effective density.Searches for decays of HNLs were performed in the l fid = 35 m long decay region (see Fig. 3) defined by the two scintillator planes SC1 and SC2, located at the distance l min = 480 m from the copper target.The decay detector covered the 3.9•10 −5 sr solid angle and had the transverse dimensions 3 × 3 m 2 , with the center displaced by 5 m from the axis.The fine-grain calorimeter at CHARM was aimed to detect inelastic scattering of electrons and muons produced in hypothetical decays of HNLs [37].The sets of tube planes P1-P5 [38] were installed to improve the reconstruction of the decay vertex and the angular resolution.
Figure 3: The layout of the CHARM facility, adopted from [19].

III. PHENOMENOLOGY OF HNLS AT CHARM A. Production
At the SPS energy of 400 GeV, HNLs with mass at the GeV scale may be produced directly either in the protontarget collisions, or in the decays of secondary particles: B, D mesons and τ leptons.The direct HNL production competes with strong interaction processes, while the production from secondary particles -with weak interactions.As a result, the latter process is dominant even taking into account small production probability of mesons [9], and the former may be completely neglected.However, similarly to the other experiment operating at SPS -NA62 in the dump mode, the CHARM experiment has no sensitivity to the HNLs produced from B mesons, implying the lower bound on the probed mass m N m Ds 2 GeV. 1et us define the HNL that mixes only with ν α by N α .
Neglecting the direct production channels, the total number of N α produced at CHARM is given by: with N cc being the total number of quark-antiquark cc pairs produced at CHARM, D i = D ± , D 0 , D s , and f c→Di the corresponding quark fragmentation fractions at SPS.The first term in the brackets describes the production from decays of D mesons (diagrams (a), (b) in Fig. 4) and the second -from τ leptons in the D s → τ → N decay chain (diagrams (c), (d) in Fig. 4).Br(D i → N α X), Br(τ → N α X) are the branching ratios.
The second term includes a small factor f c→Ds •Br(D s → τ ντ ) 5 • 10 −3 ; for the given HNL mass, it is suppressed as compared to the first term as soon as the production from D is allowed.The original analysis of the CHARM collaboration [18,19] considered the mixing α = e, µ, for which decays from D mesons are possible for any mass in the range m N < m Ds − m lα ≈ 1.9 GeV, and the production from τ decays may be completely neglected, according to the discussion above.For the τ mixing, however, the kinematic threshold of the production from D, D s → τ + N , is m Ds − m τ ≈ 190 MeV, and only the second summand in Eq. ( 3) contributes for heavier HNLs.
Let us estimate how many HNLs with τ mixing are produced as compared to those with e mixing.From (3), the ratio prod is Assuming the same values of mixing angles U 2 e = U 2 τ for the two models with pure e/τ mixing, the ratio Br(τ → N τ X)/ f c→D Br(D → N e X) varies in the 1 − 10 range for masses m N 1.3 GeV and quickly drops at the kinematic threshold m N ≈ m τ [9].In particular, for masses m N 800 MeV, where the dominant contribution to the HNL production with e mixing comes from D s , we have The mass dependence of the ratio prod obtained from Eq. ( 4) is shown in Fig. 5.It is important to note, that in the original analysis [19], as well as in the re-analyses [16,21], the production from D s has not been taken into account for the e mixing.In the mass range m N 800 MeV, this leads to the underestimate of the number of produced HNLs, N CHARM prod , by a factor 1/6 (see Fig. 4).

B. Decays and their detection
For a given number of produced HNLs, the number of detected events N (α) events for the given mixing α depends on 1. Geometrical factors -in order to be detected, produced HNLs need to point in the angular coverage of the CHARM decay volume, decay inside it, and their decay products must then reach the detector and be successfully reconstructed.These factors are: geometrical acceptance geom , i.e. the fraction of produced HNLs traveling in the direction of the CHARM detector; the mean HNL gamma factor γ N ; the decay acceptance decay , i.e. the fraction of HNL decay products that point to the CHARM detector for HNLs that decay inside the fiducial volume.

The branching ratio Br(N
ν (and their charge conjugated counterparts) used for detection at CHARM [19].The formula for N (α) events is: where is the decay probability, and det,ll is the reconstruction efficiency for the given channel.We will see below that geometrical factors are the same for e, µ and τ mixing, while the branching ratio is smaller for the τ mixing channels, as in the former case both decays via the charged and neutral currents are relevant, while in the latter only the neutral current contribute.Let us start by considering the lower bound of the sensitivity of the CHARM experiment, i.e. the minimal mix-ing angles that it may probe (the upper bound will be discussed in Sec.IV).In this regime, the decay length of the HNL cτ N is much larger than the geometric scale of the experiment, cτ • Γ(N α ), where Γ(N α ) is the total decay width, and it is convenient to rewrite Eq. ( 6) in the form where Γ(N α → l + l − ν) is the decay width into the dilepton pair ll .We will first discuss the difference in Γ(N α → l + l − ν) between the cases of e and τ mixings.Decays into dileptons occur via charged and neutral current, see Fig. 1.
For the NC mediated processes, the kinematic threshold m N > 2m e ≈ 1 MeV is mixing-independent.In contrast, for the CC mediated process for the τ mixing this threshold is m N > m τ +m e ≈ 1.77 GeV, and HNLs lighter than τ lepton may decay into dileptons only via NC.Decay widths for the processes N α → l + l − ν for m N m l + m l may be given in the unified form where the coefficients c (α) ll ν are given in Table I [9].For N e , the largest decay width is Γ(N e → µ + e − ν µ ), where only CC contributes.The width Γ(N e → e + e − ν e ) is smaller: because both NC and CC contribute in this process and interfere destructively.The smallest width is Γ(N e → µ + µ − ν e ), with the process occurring only via NC.For N τ , there is no process N τ → eµν, while in the process N τ → e + e − ν τ only NC contributes, and thus the width is smaller than for N e : Γ(N τ → e + e − ν τ )/Γ(N e → e + e − ν e ) ≈ 0.22 (10) For the decay into a dimuon pair, we have Γ( As a result, for m N m µ the ratio of the factors l,l Γ(N α → ll ν) det,ll entering Eq. ( 7) is given by l Γ(N τ → ll) det,ll l,l Γ(N e → ll ) det,ll ≈ 0.16 (11) Here and below, we use the values of the efficiencies det,ll as reported in [19] for the HNL mass m N = 1 GeV: det,ee = 0.6, det,eµ = 0.65, det,µµ = 0.75.
In the original analysis for the e mixing by the CHARM collaboration [18,19], the Dirac nature of HNLs has been assumed (the decay widths are twice smaller), and only the CC interactions have been considered.Instead of Eq. ( 11), the ratio becomes Process c Table I: The values of c ll ν in Eq. ( 8) for different decay processes.For the process Ne → e + e − νe, we also provide the value obtained if including the charged current (CC) contribution only -the assumption used in [19].
Let us now discuss geometric factors geom , γ N , decay .It turns out that they depend on the mixing pattern weakly, and as a result the geometry does not influence the relative yield of events for e and τ mixing.Indeed, as was mentioned in Sec.III A, HNLs with τ mixing are produced in decays of τ leptons, that originate from decays of D s .Since m τ m Ds , the angle-energy distribution of τ leptons is the same as of D s (and hence also other D mesons), whose decays produce HNLs with e mixing.The kinematics of the HNL production from D and τ is similar: two-body decays (a), (c) and three-body decays (b), (d) in Fig. 4 differ mainly be the replacement a neutrino or a lepton with a hadron h = π, K.However, since m h m τ,D , the replacement does not lead to the difference in the distribution of produced HNLs.In addition, heavy HNLs with masses m N 1 GeV share the same distribution as their mother particles, and any difference disappear.Therefore, the values geom , γ N for different mixing are the same with good precision.Next, HNL decays contain the same final states independently of the mixing, and decay can also be considered the same.To summarize, the ratio events is determined only by the difference in phenomenological parameters -N (α) prod and Γ(N α → ll ν): To compare with the estimate of the number of events for the e mixing made by the CHARM collaboration in [19], N CHARM events , we need to take into account their assumptions on the description of HNL production and decays (see the discussion around Eqs. ( 4) and ( 12)).The resulting ratio is IV. RESULTS Let us now derive the CHARM sensitivity to the τ mixing.In [19], it has been shown that the dilepton decay signature at CHARM is background free.Therefore, 90% CL sensitivity to each mixing is given by the condition Let us define U 2 lower,CHARM as the smallest mixing angle for which the condition ( 15) is satisfied for the assumptions of the original analysis of [19] (see the discussion above Eq.( 14)).As the number of detected events at the lower bound N (α) events scales with the mixing angle as α (where U 2 α comes from the production and another U 2 α from decay probability), we can use Eqs.( 14) and ( 4) to obtain the lower bound of the sensitivity to the τ mixing, U 2 τ,lower , by rescaling the results reported in [19]: Using the ratio N CHARM prod /N (τ ) prod from Eq. (4) (see also Fig. 5), and the ratio of decay widths from Eq. ( 12), we conclude that in the mass range m N > 200 MeV the lower bound for the τ mixing is a factor 10 − 100 weaker than the lower bound for the e mixing reported in [19].In the domain m Ds − m τ < m N < 290 MeV, we validate the rescaled bound ( 16) by comparing it with the CHARM sensitivity to the τ mixing from [20], see Appendix A. At the upper bound of the sensitivity, the dependence of the number of events on U 2 α is complicated and the sensitivity cannot be obtained by rescaling the results of [19].Therefore, we independently compute the number of decay events at CHARM for HNLs with e and τ mixing and then calculate the sensitivity numerically using Eq.(15), see Appendix A. In order to validate this estimate, we compare the resulting sensitivity for the τ mixing with the rescaled bound (16), and find that they are in very good agreement (Fig. 7).Also, we compare our estimate for the e mixing with the CHARM sensitivity to the e mixing from [19].In our estimates, we include neutral current interactions, the production from D s mesons, and assume that HNLs are Majorana particles.Due to these reasons, we find that for small mixing angles U e and above m N 1 GeV, the bound imposed by CHARM may be actually improved by up to a factor 3 − 4. Let us comment on errors of our estimates.We used the values of reconstruction efficiencies rec,ll reported in [19] for the HNL mass m N = 1 GeV.Hence, the calculation may be further refined by including HNL mass dependent reconstruction efficiencies.However, as the study [20] performed for the τ mixing and masses m N < 290 MeV has shown similar efficiency, we do not expect any significant changes.
Our final results for the τ mixing are given in Fig. 6, where we show the domain excluded by previous experiments together with updated CHARM bounds, and the sensitivity of the future experiments mentioned in Sec.I, together with SHiP [40].Comparing with Fig. 2, we find that in the mass range 290 MeV < m N < 1.6 GeV our results improve previously reported bounds on the mixing angle U 2 τ by two orders of magnitude.In particular, it excludes a large part of the parameter space that was suggested to be probed by the future experiments.For instance, Belle II, FASER, DarkQuest and the double bang signature at IceCube have sensitivity only in the narrow domain above the CHARM upper bound, while NA62 may slightly push probed angles to lower values.The same is the case for the sensitivity of KM3NeT in the regime of the atmospheric beam dump, which we estimate in Appendix B. Significant progress in testing the mixing of HNLs with ν τ can be achieved by LHCb, which probes the complementary mass range m N > 2 GeV, and dedicated intensity frontier experiments, with SHiP being optimal for searches of HNLs from decays of D mesons and τ leptons.The left panel: fraction of HNLs that point towards the detector (blue line) and fraction of HNLs whose decay products point towards the detector (red line).The middle and right panels: comparison of our estimates of the constraint from the CHARM experiment on the pure e (the middle panel) and τ mixing (the right panel), with bounds reported in [19] and [20].We show two estimates: the red line corresponds to the rescale of the bound on the e mixing from [19] (see Sec. IV for details), while the blue line is our independent estimate based on Eq. (A1).
with atmospheric particles.If having significantly large lifetimes, produced HNLs may enter the detector volume of KM3NeT located deep in the Mediterranean Sea and decay there.In order to probe the parameter space of HNLs, it is necessary to distinguish their decays from interactions of SM particles also produced in the atmosphere: neutrinos and muons.KM3NeT may only distinguish two event types: track-like, which corresponds to muons penetrating through the detector volume, and cascade-like, which originates from other particles such as electrons and hadrons.Scatterings of neutrinos inside the detector volume produce cascade-like (if no high-energy muons are produced) or combined cascade-like + track-like signature (if high-energy muons are produced), while penetrating atmospheric muons give rise to track-like signature.A possible way to distinguish these events from HNLs is to look for their decays into a di-muon pair, N → µμν τ .They produce a signature of two tracks originated from one point inside the detector volume.Detectors of KM3NeT have energy and angular resolution sufficient precise for resolving the two tracks down to energies of 10 GeV [44].Therefore, we believe that the dimuon signature may be reconstructed in the background free regime with high efficiency.

Analytic estimates: comparison with CHARM
Now, let us discuss the sensitivity of KM3NeT to HNLs.Let us first compare the amount of HNL decay events at CHARM and KM3NeT for given value of the mixing angle at the lower bound of the sensitivity using simple analytic estimates.According to Eq. ( 7), for the ratio of decay events at these experiments we have where σ pp→ccX (E p ) is the energy-dependent charm production cross-section which we use from FONLL [45] and from [46], and σ pp,total is the total pp-cross-section, which we use from [47].The integrand in (B3) is the product of two competing factors:

Figure 1 :
Figure 1: Diagrams of leptonic decays of an HNL that mixes purely with να via the charged (the left diagram) and the neutral current (the right diagram).

Figure 2 :
Figure 2: The parameter space of HNLs with the pure e (the left panel) and τ (the right panel) mixing.The current bounds are from CHARM[19,20], NA62[22], T2K[23], Belle[24,25], DELPHI[26], NOMAD[27], ArgoNeuT[28], see also[11] for a review and the references therein.For the pure τ mixing, we do not show the constraints imposed by the T2K experiment (unlike it is done in[11]), since they are reported for non-zero couplings U e/µ which dominate the production a .Constraints from the CHARM experiment are taken from the literature, while our re-analysis for them is shown in Fig.6.The light gray domain corresponds to couplings that are either excluded by BBN[12,14] or too small to provide active neutrino masses.For the pure τ mixing, we also show sensitivities of the next generation Intensity frontier experiments (see text for details).In cyan, we show HNL parameter space that may be probed by neutrino observatories: the solid line shows the sensitivity of IceCube to the "double bang" signature from[29], while the dashed line corresponds to the sensitivity of KM3NeT to decays of HNLs produced in the atmosphere, see text and Appendix B for details.

NFigure 5 :
Figure 5: The HNL mass dependence of the ratio of the numbers of produced HNLs with pure τ and e mixing N (τ ) prod /N (e) prod , see Eq. (4), assuming the same values of the mixing angles U 2 e = U 2 τ for the two models.The solid line corresponds to N (e) prod calculated keeping the production from all D mesons D + , D 0 , Ds, while the dashed line corresponds to the estimate of N (e) prod ≡ N CHARM prod calculated without the contribution of Ds, as has been done in the analysis [19] by the CHARM collaboration (see text for details).The kinks at mN = mD s −mτ ≈ 200 MeV an mN ≈ mτ −mρ ≈ 1 GeV correspond to kinematic thresholds of the production channels Ds → N + τ , τ → N + ρ correspondingly.

2 Figure 7 :
Figure7: The left panel: fraction of HNLs that point towards the detector (blue line) and fraction of HNLs whose decay products point towards the detector (red line).The middle and right panels: comparison of our estimates of the constraint from the CHARM experiment on the pure e (the middle panel) and τ mixing (the right panel), with bounds reported in[19] and[20].We show two estimates: the red line corresponds to the rescale of the bound on the e mixing from[19] (see Sec. IV for details), while the blue line is our independent estimate based on Eq. (A1).

N 2 •
Γ(N τ → ll) det,ll Γ(N τ → µµ) 10 13 (see Fig.7is the number of cc pairs detectable fraction of HNL decay events at CHARM.N KM3NeT cc is the amount of cc pairs produced in the upper hemisphere propagating to KM3NeT,N KM3NeT cc 2π × 1 km 2 × 5 years × dΦ dΩdtdSdE p • σ pp→ccX σ pp,total dE p 10 12 ,(B3) dΦ dΩdtdSdEp , which decreases with the proton's energy, and σ pp→ccX (E p ), which increases, see Fig. 8.We approximate the ratio of the mean HNL γ factors by the ratio of the mean γ factors of D mesons: where we calculate γ KM3NeT Ds using the cc distribution dΦ dΩdtdSdEp • σ pp→ccX , assuming that E D ≈ E p /2.