Search for sterile neutrinos decaying into pions at the LHC

We study the possibility to observe sterile neutrinos with masses in the range between 5 GeV and 20 GeV at the LHC, using the exclusive semileptonic modes involving pions, namely W to lepton + N to n pions + lepton+lepton (n = 1, 2, 3). The two pion and three pion modes require extrapolations of form factors to large time-like $q^2$, which we do using vector dominance models as well as light front holographic QCD, with remarkable agreement. This mass region is difficult to explore with inclusive dilepton+dijet modes or trilepton modes and impossible to explore in rare meson decays. While particle identification is a real challenge in these modes, vertex displacement due to the long living neutrino in the above mass range can greatly help reduce backgrounds. Assuming a sample of $10^9$ W bosons at the end of the LHC Run 2, these modes could discover a sterile neutrino in the above mass range or improve the current bounds on the heavy-to-light lepton mixings by an order of magnitude, $U_{l N}^2 \sim 2 \times 10^{-6}$. Moreover, by studying the equal sign and opposite sign dileptons, the Majorana or Dirac character of the sterile neutrino may be revealed.

In order to naturally explain the smallness of the observed neutrino masses, right-handed sterile neutrinos are usually introduced in the SM extensions, creating scenarios with a seesaw mechanism [8][9][10][11][12][13][14]. In the past decades, many experiments have searched for sterile neutrinos with masses in the broad range from eV to TeV [15].
Current searches for sterile neutrinos with mass m N above 100 GeV [16,17] are based on the inclusive processes pp → W * X, W * → ± ± jj [18][19][20][21]. For m N below M W , the jets tend to fall below the necessary cuts that are imposed to reduce backgrounds, so these inclusive semileptonic modes cease to be useful. Also purely leptonic modes W ( * ) → ν have been proposed [22]. The purely leptonic modes [23][24][25][26][27][28] are also an alternative to ± ± jj at the LHC, as they are not affected by the hadronic backgrounds for m N below M W , even though they have the problem of missing energy due to the undetectable final neutrino. Also, in these three-lepton modes one should not use modes with a pair of leptons with same flavor and opposite charge, such as e + e − or µ + µ − , in order to avoid the background from radiative pair production. Nevertheless, for neutrino masses below ∼ 20 GeV, the neutrino may live long enough to leave an observable displacement from its production to its decay point, a feature that will also help reduce such background [29][30][31].
In this work we consider the exclusive semileptonic processes W → N , N → π, ππ and πππ, which are modes with no missing energy. Among these modes, we do not consider N → e ± π ∓ in the secondary process to avoid misidentification of electrons and pions. The most promising modes should be W + → µ + N , N → µ + nπ (n = 1, 2, . . .) for a Majorana sterile neutrino, and W + → e + N , N → µ − nπ for a Dirac neutrino.
The exclusive single-pion mode N → π is suppressed with respect to the inclusive process, modeled by the production of an open quark pair N → qq , by a factor f 2 π /m 2 N , where f π is the pion decay constant. This suppression is clearly milder for lower masses of the sterile neutrino. Moreover, for lower m N the sterile neutrino width is also smaller, giving time for the neutrino to travel a measurable distance before decaying. This effect will show as a spatial vertex separation in the detector, a distinctive feature that helps reduce backgrounds [29,30].
Less suppressed than the decay into a single pion is the decay into two pions, i.e. N → ππ, or even into three pions, because the two quarks from the weak vertex are not forced to be bound into a single pion wave function [32]. Although there is a suppression from the 3-particle or 4-particle phase space, it is not as strong as the wave function restriction of the single pion case.
As stated above, at the LHC the inclusive mode ± ± jj is sensitive to neutrino masses m N above 100 GeV and the purely leptonic ± ± ∓ ν modes cover the region below M W down to about 20 GeV. Complementary, for m N below 5 GeV, rare decays of B, D and K mesons or tau leptons are more appropriate [33][34][35][36][37]. In this work, we focus on neutrino masses within the window from 5 GeV to 20 GeV, using the exclusive processes W + → + N followed by the pionic decays N → π − + , N → π 0 π − + and N → π − π − π + + . These processes can give the LHC the ability to explore the corresponding mass region 1 .
This article is organized as follows. In section. II, we calculate the pionic decays of a sterile neutrino. In section. III, we detail our numerical results and discussions. We conclude in section. IV.

II. PIONIC DECAYS OF A STERILE NEUTRINO
The leptonic sector in a generic SM extension includes one or more extra neutral lepton singlets, N , in addition to the three generations of left-handed SM SU (2) L lepton doublets.
The neutral lepton singlets N are known as sterile neutrinos, because they do not directly interact with other SM particles in the absence of any mixing with the active neutrino sector.
After electroweak symmetry breaking, the gauge interaction of the sterile neutrino with the SM fields can be written as where g denotes the weak coupling constant and U N the mixing matrix between the active and sterile neutrinos. Here the sterile neutrino can be a Dirac or Majorana fermion.
At the LHC, sterile neutrinos with masses around 5 ∼ 20 GeV will be mainly produced from the decay of on-shell W bosons. From the lagrangian in Eq. (1), the decay rate W + → + N can be easily calculated. Neglecting the lepton mass, the branching ratio B(W + → + N ) is: where Γ W 2.085 GeV is the total decay width of the W boson [38]. From here, the heavy neutrino N can decay in several modes, depending on its mass. Here we are interested in the decays into pions, namely N → π ∓ ± , N → π 0 π ∓ ± and N → π ∓ π ∓ π ± ± . Both charged modes will occur for a Majorana N , while for a Dirac N only the N decays into a negative charged lepton will be produced.
A. The decay N → π − + : The two-body decays of a heavy neutrino have been studied in detail [20]. The mode N → π − + is a charged current process: For an unpolarized N , the decay into the opposite charge mode N → π + − is the same as above.
In an obvious notation, m π − and m N denote the mass of the charged pion and sterile neutrino, respectively; V ud is the CKM matrix and f π is the pion decay constant; the function λ(x, y, z) is defined as λ(x, y, z) = x 2 + y 2 + z 2 − 2(xy + yz + zx). Now, the formation of a single pion in the final state is relatively suppressed with respect to multi pion modes, because it requires the two produced quarks to remain close together. Indeed, we can compare the above expression with that of open quark production, which is an estimate of the inclusive rate. Neglecting the masses in the final state, the rate is Γ(N → qq ) ∼ (G 2 F /64π 3 )|V qq | 2 |U N | 2 m 5 N , and therefore the suppression factor of the exclusive π mode is approximately: For m N = 10 GeV, this suppression is R ∼ 0.6%. Consequently, in order to have a larger rate, we should also consider the decays with two and three pions in the final state.
B. The decay N → π 0 π − + : The decay into two pions, N → π 0 π − + , is similar to the tau lepton decay τ − → π 0 π − ν τ in terms of their interaction lagrangian and Feynman diagram, except for the lepton flavor and charge. Indeed, the hadronic current in both decays has the same expression in terms of form factors. However, one must be aware that the kinematic range for the form factor in the N decays is extended to higher q 2 , so an extrapolation of the tau form factor will be required.
Considering the above, the differential decay rate for N → π 0 π − + can be written as The decay rate is then obtained after integrating over s: within the limits s − = (m π − + m π 0 ) 2 and s + = (m N − m ) 2 . It is easy to check that, after the replacement (N, + ) → (τ − , ν τ ), our expression coincides with that of the decay width for [39,40]. The main challenge here is to estimate the form factor F − (s). For low s, the form factor can been calculated in the framework of Chiral Perturbation Theory (ChPT) [41]. In the time-like region, i.e. s > 0, F − (s) is experimentally known from τ − → π − π 0 ν τ decay [42], but only in the limited range 2m π < √ s < m τ . On the other hand, the neutral form factor F 0 (s), defined by which is related to F − (s) by conservation of the vector SU (2) isospin current (CVC) as F 0 (s) = F − (s), can be obtained from current e + e − → π + π − data [43] up to √ s < 3 GeV. In our case, we need to know the form factor up to m N ∼ 20 GeV, so we will need to do an extrapolation.
For this purpose we consider two theoretical models: one is the vector meson dominance model (VDM) used by the BaBar collaboration [43], and another is based on light front holographic QCD (LFH) [44]. The VDM parametrization for F 0 (s), which includes ρ-ω mixing, is given by [43] F 0 (s) = 1 where the Breit-Wigner functions BW i (s, m i , Γ i ) and numerical values of all the coefficients fitted to the e + e − → π + π − data in the region √ s < 3 GeV, can be found in [43] (see also [45]).
Here no vector mesons heavier than 3 GeV contribute to F 0 (s), because no light unflavored vector mesons heavier than 3 GeV have been observed yet, qc and qb mesons (e.g. D * (2007) and B * (5324)) do not decay to ππ through strong interactions and are therefore suppressed, and the decays of cc and bb mesons (e.g. J/ψ and Υ) into two pions are OZI suppressed.
Since isospin is broken by the small quark mass difference m d − m u , the mesons ρ 0 (I = 1) and ω (I = 0) are not exact isospin eigenstates, but they mix with each other [46]. This feature results in a small structure near √ s = 0.78 GeV in the e + e − → π + π − spectrum. For the form factor F 0 (s), this effect is included in Eq. (9). However, since there is no charged ω meson, no ρ − ω mixing should be present in F − (s) and, moreover, since we are obtaining F − (s) from the CVC isospin relation F − (s) = F 0 (s), isospin violation effects should be omitted. Nevertheless, since the parameters in F 0 (s) were found by fitting the data altogether and this isospin violation effect is very small, we kept the ω term in F − (s).
Concerning the holographic QCD formalism [44], the form factor is dominated by the twist-2 and twist-4 terms (in this formalism, the twist 'τ ' corresponds to the number of partons in the light-front Fock state and is also related to the number of poles in the form factor, which is τ − 2). In this approximation, the form factor is: where γ = 12.5 % is the probability of the twist-4 component. The two components are: and and the following parameters were used for their best fit up to s ∼ 3 GeV 2 : m ρ = 0.775 GeV, According to Ref. [44], their fit agrees with data up to s ∼ 6 GeV 2 . The pion form factor according to these two models is shown in Fig. 1 in both models one finds the QCD behavior sF (s) → constant at large s, as shown in Fig. 1 (right). Since our purpose is to estimate N → ππ for m N within the range 5 GeV to 20 GeV, it is necessary to have a good estimate of the form factor extrapolated beyond current data.
Fortunately the extrapolated region (large s) is where theory is more reliable, and the most difficult part to address theoretically is the resonance region where data is available. Indeed, as the above study shows, the extrapolation with two different treatments give quite similar results, so we expect this estimate to be reliable for the purpose of our work.
In much the same way as in the two-pion mode of Section II.B, the differential decay rate of the general hadronic decay N → h 1 h 2 h 3 + can be written in terms of form factors with an expression identical to that of the tau decay τ − → h 1 h 2 h 3 ν τ [47], again provided that the form factors are extrapolated to larger values of q 2 . Let us first denote the momentum and mass of the hadron h i (i = 1, 2, 3) by p i and m i respectively, and define the momentum of the hadronic part q µ = (p 1 + p 2 + p 3 ) µ . The differential decay rate is then: Here the functions ω i (q 2 ), i =A, B, SA, are defined as [48]: with the variables s 3 = (p 1 + p 2 ) 2 and s 2 = (p 1 + p 3 ) 2 , and the integration limits: The integrands W i (q 2 , s 2 , s 3 ) are the following Lorentz invariants: Here V µ 1 , V µ 2 and V µ 3 are the kinematical vectors: while F A 1 (q 2 , s 2 , s 3 ), F A 2 (q 2 , s 2 , s 3 ), F V 3 (q 2 , s 2 , s 3 ) and F P 4 (q 2 , s 2 , s 3 ), are the form factors that appear in the hadronic current matrix element for the decay N → h 1 h 2 h 3 + [47]: In this parametrization, the form factors F A 1 and F A 2 correspond to axial vector transitions (J P = 1 + ), while F V 3 and F P 4 to vector (J P = 1 − ) and pseudoscalar (J P = 0 − ) transitions, respectively. Now, in the N → π − π − π + + decay, the partial conservation of the axial current (PCAC) implies that the pseudoscalar form factor F P 4 is proportional to m 2 π /q 2 , and therefore suppressed in most of the kinematical region [49][50][51]. Moreover, in the isospin symmetry limit, G-parity conservation leads to F V 3 = 0 [48]. Therefore, we only need to consider the two axial form factors F A 1 and F A 2 . For q 2 < m 2 τ these two form factors were first obtained from the CLEO τ decay data, and then can be extrapolated to the q 2 > m 2 τ region by using an appropriate parametrization. There are many parametrizations of the form factors in the literature. We chose the so-called CLEO parametrization [52], which is used in the TAUOLA package [53,54]. This parametrization includes the transitions between various intermediate resonances a 1 (1260)/a 1 (1640) → π + f 0 (500), f 2 (1270), f 0 (1370), ρ(770), and ρ (1450). Just as in the N → π 0 π − + decay, meson resonances heavier than the τ lepton do not contribute to the form factors. At s 2 = s 3 = 1 GeV 2 , the numerical results of F A 1 (q 2 , s 2 , s 3 ) and F A 2 (q 2 , s 2 , s 3 ) are shown in Fig. 2.
The decay rate is obtained after integration over q 2 : where the integration limits are: One can check that, after the replacement (N, + ) → (τ − , ν τ ), our formulae coincide with those for the tau lepton decay rate τ − → h 1 h 2 h 3 ν τ given in [47].

III. NUMERICAL ANALYSIS
With the theoretical framework described in the previous section, we want to estimate the exclusive semileptonic decay rates of N into π , 2π and 3π , for a neutrino N of mass in the range 5 to 20 GeV, produced at the LHC in the process W → N . These exclusive modes could be observed in pp collisions provided the pions can be identified and the background can be reduced considering the spatial displacement between the production and decay vertices of N . This displacement could be observable for m N below 20 GeV [29,30]. We first study the invariant mass distributions for the two-pion and three-pion modes, for two benchmark neutrino masses, m N = 5 GeV and 15 GeV, and then study the integrated rates of the single pion, two-pion and three-pion modes as a function of m N . Within the quark-hadron duality [55], the open quark process N →ūd + is considered to be the same as the inclusive semileptonic process. As shown in Fig. 3, the exclusive mode is much smaller than the inclusive process in the invariant mass region far from the ρ meson peak, while they are of comparable magnitude in the region near the peak. Moreover, considering the integration over √ s in these figures, one can notice that the total decay width of N → π 0 π − + is predicted, as it should, to be smaller than the inclusive estimate N →ūd + , and this difference is enhanced for larger m N , as there will be more modes with larger hadron multiplicities as N is heavier. . Canonical decay rates for N → π − + , N → π 0 π − + , N → π − π − π + + and N →ūd + (left) and canonical branching ratios for the full processes W + → + + nπ (right) as a function of the neutrino mass m N , with all mixing factors |U N | removed. To obtain the actual values, the canonical values must be multiplied by the factor |U N | 2 (left), or |U N | 4 / l |U N | 2 (right).
In Fig. 4 (left) we show the "canonical" decay rates for the modes N → + nπ and the inclusive estimate given by N →ūd + as a function of the neutrino mass m N ("canonical" here means that all lepton mixing elements |U N | are removed from the expressions). In Fig. 4 (right) we show the "canonical" branching ratios for the full processes W + → + N → + + nπ and W + → + N → + +ū d. The actual rates (left) and branching ratios (right) can be obtained by multiplying these canonical values by |U N | 2 (left) and |U N | 4 / l |U N | 2 (right), respectively.
These factors can be easily deduced considering the expressions for the partial N decay rates in Eqs. (3), (5) and (13), and the expression for the neutrino width Γ N [20,23,56]: which is a good approximation in the range 5 GeV< m N < 20 GeV.
From Fig. 4, it can be seen that Γ(N → 3π )/Γ(N → 2π ) ≈ 0.7 . This is expected, as modes with more particles in the final state tend to be smaller than those with fewer particles. However, Fig. 4 also shows that the single pion mode is suppressed with respect to the two-pion mode, Γ(N → π )/Γ(N → 2π ) ≈ 0.4, thus contradicting the previous argument, but, as explained in the Introduction, the probability of producing a single pion from a pair of energetic quarks is small and constitutes an even stronger suppression that the two-or three-pion final state.
We would like to remark that the quark level process N →ūd + can be considered quite confidently as the inclusive semileptonic process -what is usually called the quark-hadron duality [55]-for m N > 5 GeV, because the pion multiplicity can be quite large for these invariant masses, and the region is well above the resonance region. This is not the case in b → c transitions at invariant masses near 5 GeV, where the hadronic modes are very few.
Indeed, the recent BaBar and Belle data indicate that a sizable duality violation occurs in the [57]. Finally, let us estimate the expected number of W → N → nπ events. According to Ref. [58], at the end of the LHC Run II one may expect a sample of N W ∼ 10 9 W decays. Considering the canonical branching ratios in and B(W + → π + π − π − µ + µ + ) ≈ 0.7 B(W + → π 0 π − µ + µ + ). Including all the modes, we have In order to obtain more than 5 events, we must have: GeV< m N < 50 GeV come from DELPHI [59], which are |U eN | 2 , |U µN | 2 , |U τ N | 2 2.1 × 10 −5 .
Consequently, using the modes W ± → nπ ± ± , a sample of 10 9 produced W bosons may be able to improve the current upper bound on |U N | 2 by nearly an order of magnitude. Moreover, these pionic modes seem to be amongst the most promising to explore the sterile neutrino mass The modes with opposite charge leptons, i.e. e + e − , µ + µ − , µ + e − etc. will also occur and with similar rates as those above, provided N is Majorana. Although the opposite charge dilepton modes suffer from other backgrounds, they may help improve the statistics for a Majorana N .
On the other hand, if N were a Dirac neutrino, all the equal-sign dilepton modes we analyzed will be absent, but the opposite charge dileptons will still be present and with the same rates as calculated above.
We find, as expected, that the single pion mode is suppressed with respect to the two-pion and three-pion modes, while the latter two rates are comparable. This suppression is explained as the low probability of forming a single pion from the two energetic light quarks produced in the decay of N with mass of a few GeV or more, and expressed as the small fraction f π /m N .
In the two-pion mode, the calculation involves an extrapolation of the pion form factor at large time-like q 2 . For this purpose we use two models of form factors, one based on the vector meson dominance model (VDM) used by the BaBar collaboration and another based on light front holographic QCD (LFH). The two treatments show a remarkable agreement in the extrapolated region, which may be considered both an indication of the success of these two treatments as well as a confidence on our estimate.
In the three-pion mode, we use an extrapolation from tau lepton decays based on the socalled CLEO parametrization, which includes a series of resonances. This mode gives a value that is consistently lower than the two-pion mode by a factor of ∼ 0.7 for all the m N mass range, which is also a sign of confidence on the model estimate.
Assuming ∼ 10 9 W's produced by the end of the run II of the LHC, and in the bestcase scenario of negligible background, our studied pionic modes could be able to put limits to |U N | 2 ∼ 2 × 10 −6 , one order of magnitude more stringent than the bounds coming from DELPHI. Moreover, the equal sign dilepton modes will provide signals or bounds for Majorana N , while the study of opposite sign dileptons will provide bounds for a Dirac N , provided the equal sign modes are not seen. As a final comment, while we did not study the backgrounds, which may turn to be very difficult to eliminate considering the pion identification, one should still keep in mind that for sterile neutrino masses below ∼ 20 GeV, the neutrino may live long enough to cause an observable displacement between the primary and secondary lepton, a feature that may greatly help reduce the backgrounds.
After finishing this study we learned of two recent works [60,61] that deal with the issue of displaced vertices in N related processes.