Search for heavy neutral leptons in $W^+\to\mu^{+}\mu^{\pm}\text{jet}$ decays

A search is performed for heavy neutrinos in the decay of a $W$ boson into two muons and a jet. The data set corresponds to an integrated luminosity of approximately $3.0 \text{ fb}^{-1}$ of proton-proton collision data at centre-of-mass energies of 7 and $8 \text{ TeV}$ collected with the LHCb experiment. Both same-sign and opposite-sign muons in the final state are considered. Data are found to be consistent with the expected background. Upper limits on the coupling of a heavy neutrino with the Standard Model neutrino are set at $95\%$ confidence level in the heavy-neutrino mass range from 5 to $50 \text{ GeV}/c^2$. These are of the order of $10^{-3}$ for lepton-number-conserving decays and of the order of $10^{-4}$ for lepton-number-violating heavy-neutrino decays.


Introduction
Many theories beyond the Standard Model (SM) predict the existence of heavy neutral leptons (HNLs) to explain the smallness of neutrino masses [1][2][3]. These leptons, N , could be observed at collider experiments if their masses are at the electroweak scale. The HNLs may mix with the light SM neutrinos ν , with a strength governed by the coupling V N . The mixing matrix is not expected to be flavour diagonal, which leads to signatures with transitions between different lepton flavours. Experimentally, direct searches for a generic heavy neutrino are performed through their mixing with each flavour of SM neutrino. The HNL is expected to be long-lived if the coupling is small enough. This analysis searches for the mixing of a heavy neutrino with a muon neutrino, taking advantage of the high reconstruction efficiency for muons at LHCb. The HNL mass range covered is between 5 and 50 GeV/c 2 . The dominant HNL production mechanism in this mass range is via the decay of gauge bosons, W ± → ± N and Z → νN . Both lepton-number-violating and lepton-number-conserving decays of a heavy neutrino are considered. The heavy neutrino is assumed to have negligibly small lifetime.
The DELPHI collaboration was first to set a limit on these types of decays considering Z → νN decays in e + e − collisions at the Z resonance, where both long-and short-lived signatures were analysed [4]. The upper limit on the Z → νN branching fraction of 1.3 × 10 −6 at 95% confidence level (CL) for N masses between 3.5 and 50 GeV/c 2 leads to one of the most stringent constraints on the coupling in this mass range. At the LHC, a more promising detection approach for N with mass below the weak scale are leptonic decays of W bosons, W ± → ± N . Searches by the ATLAS [5-9] and CMS [10][11][12][13][14][15][16] collaborations at centre-of-mass energies of 7, 8 and 13 TeV typically probed larger neutrino masses, from 40 GeV/c 2 up to 2700 GeV/c 2 , employing a signature of two same-sign leptons and two jets. A recent search performed by the CMS collaboration also included final states with at least one jet, extending the probed heavy-neutrino mass range down to 20 GeV/c 2 [17]. In the mass range studied in this analysis, searches of promptly decaying heavy neutrinos in leptonic final states of the W boson at centre-of-mass energy of 13 TeV by the ATLAS [18] and CMS [19] collaborations set constraints comparable to that of the DELPHI collaboration. A long-lived signature has also been explored by the ATLAS collaboration, excluding coupling strengths down to about 10 −6 between 4.5 and 10 GeV/c 2 , and hence representing the most stringent limit to date in this mass range [18].
The branching fraction (B) for the decay of a W boson into a muon and a heavy neutrino is suppressed with respect to the SM decay W + → µ + ν by the mixing of the active neutrino with the heavy neutrino and a phase-space factor according to Ref. [20] The heavy neutrino decays via neutral or charged current interactions N → νZ ( * ) or N → µ ± W ∓( * ) , where the Z and W bosons can be on-or off-shell. The corresponding branching fractions are computed based on Refs. [21,22]. The total width is given by the sum of the partial decay widths of charged and neutral current interactions. If the neutrino is a Majorana particle, an additional lepton-number-violating decay contributes to the same final state, with the same partial decay width as the lepton-number-conserving decay. The branching fraction to any non-charge-specific final state is unaffected, but the lifetime is a factor of two smaller than if the neutrino were a Dirac particle.   [21,22]: (left) the branching fractions to final states with a muon and (right) the lifetime, assuming a coupling of 10 −5 (current limit).
The left plot of Fig. 1 shows the branching fraction for HNL decay modes with a muon in the final state as a function of the heavy-neutrino mass. The difference between the HNL decay modes to quarks is mainly due to CKM matrix elements [23,24], with the quark masses only playing a minor role at low heavy-neutrino masses. The branching fraction of the decay N → µµν is about one order of magnitude smaller than that of the N → µqq mode, due to negative interference between charged and neutral current interactions. In the right plot of Fig. 1 the lifetime is shown as a function of the heavy-neutrino mass assuming a coupling of 10 −5 , which is the current bound in the mass range of interest. In the low-mass regime, the lifetime is of the order of 10 ps, while at higher masses the lifetime is so small that the decay can be considered prompt.
In this paper, a search is presented for a prompt HNL in the decay 1 W + → µ + N with N → µ ± qq , as depicted in Fig. 2. Data collected by the LHCb experiment in proton-proton collisions at centre-of-mass energies of 7 TeV in 2011 and 8 TeV in 2012 are used, corresponding to integrated luminosities of 1.0 and 1.9 fb −1 [25], respectively.
The experimental signature consists of two muons and one or two jets depending on the HNL mass. The muon from W decay, denoted as µ W , carries significant transverse momentum, while the muon from N decay, denoted as µ N , has lower momentum. Both same-sign and opposite-sign muons are considered, allowing for the possibility that the HNL has a Majorana nature. The signal yields for both categories and several mass hypotheses in the range 5 − 50 GeV/c 2 are extracted from the data and normalized with respect to the W + → µ + ν decay. Corresponding upper limits are then set on the product of coupling and branching ratio.
The paper is organised as follows. In Section 2 the detector, data and simulation samples are described, and in Section 3 the selection of signal and normalisation candidates is discussed. Section 4 contains the results and conclusions are drawn in Section 5.

Detector and simulation
The LHCb detector [26,27] is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing b or c quarks. The detector includes a high-precision tracking system consisting of a siliconstrip vertex detector surrounding the pp interaction region [28], a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about 4 Tm, and three stations of silicon-strip detectors and straw drift tubes [29] placed downstream of the magnet. The tracking system provides a measurement of the momentum, p, of charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV/c. The minimum distance of a track to a primary proton-proton collision vertex (PV), the impact parameter (IP), is measured with a resolution of (15 + 29/p T ) µm, where p T is the component of the momentum transverse to the beam, in GeV/c. Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov detectors [30]. Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillating-pad (SPD) and preshower detectors, an electromagnetic and a hadronic calorimeter. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers [31]. The online event selection is performed by a trigger [32], which consists of a hardware stage, based on information from the calorimeter and muon systems, followed by a software stage, which applies a full event reconstruction. For the events selected for this analysis, the trigger requires at least a single muon with p T > 1.48 (1.76) GeV/c at the hardware stage in 2011 (2012), and includes an upper threshold of 600 hits in the SPD to prevent high-particle multiplicity events from dominating the processing time. A muon with p T > 10 GeV/c is required at the software stage.
Simulated samples were generated for the signal decay with both opposite-and same-sign muons in the final state, in equal amount. Samples were generated for HNL masses of 5, 10, 15, 20, 30, and 50 GeV/c 2 , using the minimal mixing scenario model [33]. The parton level process is generated with MadGraph 5 [34,35], while Pythia 8 [36], with a specific LHCb configuration [37], is used for the generation of the underlying event, fragmentation and hadronisation. Decays of hadronic particles are described by EvtGen [38], in which final-state radiation is generated using Photos [39]. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [40] as described in Ref. [41]. Simulated background samples are generated using Pythia 8. The DYTurbo [42] program is used to correct the kinematic distributions of the simulated W + → µ + ν background.

Event selection
Signal candidates are reconstructed from a pair of charged tracks identified as muons and a single jet. First, the high-momentum muon, µ W , is selected. Both the hardware and software trigger decisions are required to be associated to the high-momentum muon candidate. The track is required to have transverse momentum larger than 20 GeV/c, to be of good quality and to have a high significance of the track curvature to remove high transverse momentum tracks with poorly determined charge. The high-momentum muon candidate is also required to have small relative energy deposition in the calorimeters to reject pions and kaons misidentified as muons. The muon selection criteria for the normalisation channel W + → µ + ν are the same as for the high-momentum muon of the signal.
The lower-momentum muon candidate, µ N , is required to have transverse momentum higher than 3 GeV/c. The combined invariant mass of the µ N and µ W candidates must be in the range 20 to 70 GeV/c 2 to suppress the background from Z → µµ decays. Depending on the relative charge of the two muons the candidates are classified as same-sign (SS) or opposite-sign (OS).
Jets are reconstructed following a particle flow approach [43], using tracks of charged particles and calorimeter energy deposits as inputs. To prevent overlap between jets and signal muons, tracks identified as muons with a transverse momentum greater than 2 GeV/c are excluded from the jet reconstruction. The anti-k T jet clustering algorithm is used [44], with a distance parameter R = (∆φ) 2 + (∆η) 2 = 0.5, where φ is the azimuthal angle and θ the pseudorapidity. The jet four-momentum is calculated from the four-vectors of its constituents, and corrected for pollution from pile-up and the underlying event using the per-event particle multiplicity [43]. To enhance the jet purity the fraction of the jet energy carried by charged particles should be at least 0.1, the jet must have p T > 10 GeV/c and contain at least one track with p T > 1.2 GeV/c. Only candidates with at least one jet passing these criteria are retained. Jets are combined with lower-momentum muon candidates to form N → µ N jet candidates, which are required to have invariant mass smaller than 80 GeV/c 2 and a transverse momentum greater than 10 GeV/c. The selected heavy-neutrino candidates are then combined with a high-momentum muon candidate to form W candidates. Since the assignment of the two muons is ambiguous if they both satisfy the high-momentum muon selection, the mass, m(µ N jet), of the µ N jet combination is required to be smaller than that of the µ W jet combination. Only the candidates within 20 GeV/c 2 of the known W mass [45] are retained.
A scale factor is applied to the jet four-momentum, constraining the invariant mass of the µ W µ N jet system to the known mass of the W boson. This leads to a significant improvement in the resolution of m(µ N jet) and diminishes the sensitivity of the heavyneutrino mass distribution to the jet energy scale. Dominant background sources are charged weak currents, in particular pp → W + X with W → µν or W → τ ν, neutral electroweak Drell-Yan processes pp → γ/Z ( * ) + X with γ/Z ( * ) → µµ, τ τ , heavy flavour bb → Xµ and cc → Xµ, and Xµ production from light QCD (u, d, s). In the same-sign muon channel the Drell-Yan type background contributions are highly suppressed, while in the opposite-sign muon channel the contribution from low-mass Drell-Yan processes remains a prominent irreducible background.
The background is suppressed by placing requirements on the IP of both muons and by using three different multivariate classifiers based on a Boosted Decision Tree (BDT) algorithm [46][47][48]. The three classifiers are referred to as the µ W uBDT, the µ N uBDT and the kinematics uBDT: the first two classifiers are dedicated to the identification of the respective muons, while the latter exploits the event kinematics to distinguish the signal from the remaining background. All three are trained minimising the dependence of the signal efficiency on the true neutrino mass using the uBoost method [49]. The training of all classifiers uses a cross-validation technique [50]. The classifiers are trained using simulated decays of the heavy neutrino with same-sign muons in the final state as a proxy for signal. Both charged and neutral weak background contributions have a muon in the final state with similar kinematics to the signal high-momentum muon. The µ W classifier discriminates between the signal and heavy flavour background. It is trained using data candidates where both muons have large impact parameters (IP(µ W ) > 40 µm, IP(µ N ) > 100 µm) as a proxy for background. For both the µ N uBDT and the kinematics uBDT, a combination of the dominant background sources from simulation is used. The input variables used for each of the muon identification classifiers are the combined particle identification information from the RICH, calorimeter and muon systems, the ratio of the energy deposited in both calorimeters to the measured track momentum, and observables describing the isolation of the tracks. Additional isolation variables of different cone sizes are included among the inputs for the µ N uBDT classifier. The input variables of the kinematics classifier comprise the angular distance R between the µ N and the jet, the angle between the two muons in the rest frame of the heavy neutrino, the transverse component of the sum of the four-momentum of all particles used as particle flow input, the dimuon mass, the combined invariant mass of the dimuon and the jet, and the jet transverse momentum. The optimal requirement on the output of each BDT classifier is selected by maximising the Punzi figure-of-merit [51] for three units of significance. This is first evaluated for the µ W uBDT, followed by the simultaneous optimisation of the µ N and kinematics uBDTs. The optimal requirements are found to be the same for all the simulated signal samples. The selection is optimised for the same-sign muon signal, but it is verified to be optimal for the opposite-sign category as well.
The background sources are studied and evaluated in three control regions: one enhanced in electroweak W background components, one in heavy flavour background components and one in light QCD background components, indicated as W , bb and QCD regions, respectively. The requirements defining the control regions with respect to the signal region are reported in Table 1. An additional region, denoted as the Z → µµ region, is defined by the following criteria: both muons are required to have transverse momentum greater than 20 GeV/c and IP smaller than 40 µm and the invariant mass of the muon pair must be between 60 and 120 GeV/c 2 . In each control region the predicted

Fit strategy and results
The product of the branching fraction B(N → µ jet) and the squared coupling |V µN | 2 is proportional to the number of signal candidates, N sig , and can be written with respect to the number of W → µν candidates as where N sig and N norm denote the yields of the signal and normalisation channels and ε sig and ε norm their efficiencies. The phase-space suppression factor and the coupling term arise from the heavy-neutrino production process described by Eq. 1. The W + → µ + ν branching fraction in Eq. 1 cancels with the normalisation channel.

Normalisation channel
The yield of the normalisation channel is determined using a binned maximum-likelihood fit to the muon transverse momentum distribution separately for each year of data taking and in eight bins of muon pseudorapidity. The fit is performed separately for positively and negatively charged muons to account for the difference in production rate at LHCb. The main background contributions are γ/Z ( * ) → µµ decays and hadron misidentification (denoted as QCD). Minor contamination from Z → τ τ , W → τ ν and bb processes is also present. The templates are obtained in bins of pseudorapidity from simulation for each component, with the exception of the QCD background templates that are determined from a control sample characterised by large energy deposits in the calorimeters. The yields for the minor background contributions are fixed to their expected values from simulation. The yield for the Z → µµ component is constrained to the value obtained from the corresponding control region extrapolated according to simulation. The distribution of the muon transverse momentum for 2012 data integrated over pseudorapidity is shown in Fig. 3 with the fit result superimposed. Systematic uncertainties on the normalisation yield are estimated separately for the 2011 and 2012 data sets by varying the shape and normalisation of the templates. Replacing the QCD template with an exponential distribution and varying the W → µν templates each yield a difference with respect to the default fit of about 1%, which is assigned as a systematic uncertainty. The ratio of measured QCD yields per pseudorapidity bin between  positively and negatively charged muons deviates from unity. A systematic uncertainty of 0.7% is assigned to account for the difference with respect to the default fit when the normalisation of the QCD component is fixed bin by bin to the average of the yields. Systematic uncertainties of less than 0.1% are assigned for each of the components whose yield is fixed in the fit to account for the largest variation observed with respect to the default fit when each yield is changed by one standard deviation. For the 2011 data set an additional source of uncertainty is considered to account for the difference in templates between 2011 and 2012 simulation, resulting in 0.7% assigned systematic uncertainty. The total systematic uncertainty on the normalisation yield is 1.8% and 1.6% for the 2011 and 2012 data sets, respectively. The total yield for the normalisation channel W → µν is (795 ± 1 ± 15) × 10 3 for 2011 data and (1719 ± 2 ± 27) × 10 3 for 2012 data, where the first uncertainty is statistical and the second systematic. The total yield comprises 57% W + decays and 43% W − decays. The ratio of the measured yields for positively and negatively charged muons as a function of pseudorapidity is in good agreement with the simulation and the measurement of Ref. [52].

Efficiency ratio
The efficiency and corresponding uncertainties of the selection requirements for both the normalisation and signal samples are determined separately for each year of data taking using simulation. Corrections to account for mismodelling in simulation are derived from control samples, such as Z → µ + µ − and J/ψ → µ + µ − , and are applied to the efficiencies related to the reconstruction of the two muons, the required number of hits in the SPD and the µ W uBDT and µ N uBDT criteria. When sufficient data is available, the corrections are evaluated in bins of pseudorapidity and momentum or transverse momentum of the muon. The dominant source of systematic uncertainties arises from the different detector response 24 ± 1 ± 2 25 ± 1 ± 3 22 ± 1 ± 2 21 ± 1 ± 1 10 24 ± 1 ± 2 24 ± 1 ± 2 21 ± 1 ± 2 19 ± 1 ± 2 15 25 ± 1 ± 3 26 ± 1 ± 3 24 ± 1 ± 2 23 ± 1 ± 2 20 29 ± 1 ± 4 28 ± 1 ± 3 26 ± 1 ± 4 25 ± 1 ± 3 30 32 ± 1 ± 3 32 ± 1 ± 4 29 ± 1 ± 4 30 ± 1 ± 3 50 61 ± 3 ± 3 55 ± 2 ± 4 43 ± 2 ± 4 43 ± 2 ± 5 to jets between simulation and data. The energy scale is modelled to an accuracy of about 5%, driven mainly by the response to neutral particles, while the jet energy resolution is modelled in simulation to an accuracy of about 10% [43,53,54]. The corresponding systematic uncertainties on the efficiency ratio are evaluated in simulation by varying the jet energy by 5% for the former and by smearing the jet transverse momentum by 10% for the latter. Both resulting uncertainties vary between 5% and 11% depending on the heavy-neutrino mass, where the fluctuation is due to the limited size of the simulated samples. The overall uncertainty due to jet identification requirements, which amounts to 1.7%, is taken from Ref. [55]. A systematic uncertainty to account for the mismatch between simulation and data of the missing transverse momentum in the event varies between 1% and 2.5% depending on the heavy-neutrino mass. The uncertainties related to the efficiency of the µ W selection largely cancel for the signal and normalisation modes, since their selections are identical. The relative uncertainty on the correction factors is of the order of 2%. The ratios of efficiencies between the normalisation and signal channel, for different heavy neutrino masses, are reported in Table 2.

Neutrino mass model
The signal yield for each heavy-neutrino mass hypothesis is determined from a binned maximum-likelihood fit to the invariant mass m(µ N jet). In the fits, the normalisation channel yield, the efficiency ratio, and background yields are Gaussian-constrained to their expected values within uncertainties. The yields for the main background components are determined in the respective control regions. The yields for W → µν and Z → µµ background contributions are obtained from a binned maximum-likelihood fit of the invariant mass m(µ N jet) in the W region, and for the bb background in the bb region. The fits in the control regions are performed separately for positively and negatively charged µ W and per year of data taking with templates obtained from simulation. The expected background yields in the signal region are determined by scaling the fitted yields according to simulation. The light QCD contribution in the signal region is estimated with a different method. The efficiency of the µ W uBDT requirement ε QCD is evaluated using the normalisation channel, assuming that it factorises from the other selection criteria that suppress the QCD background. The number of light QCD events in the QCD region is obtained by subtracting from the total number of events the expected yields for the W , Z and heavy flavour background components. The result is scaled by the ratio ε QCD /(1 − ε QCD ) to determine the number of light QCD events in the signal region. The estimated background yields in the signal region are collected in Table 3 for same-sign and opposite-sign muons in Run 1 (2011 and 2012 combined) data. The uncertainty is dominated by the limited size of the simulated samples. The background predictions are tested in validation regions. These are defined by inverting one by one the requirements of Table 1 defining the signal region. The results are found to be in good agreement with the data. The templates for both signal and background contributions are determined from simulation. The light QCD background is assumed to have the same shape as the bb background and therefore a single component is included in the fit for both.

Results
The number of events observed in data in the signal region amounts to 8 and 2147 for samesign and opposite-sign muons, respectively. A single fit to the Run 1 data is performed since the 2011 and 2012 templates are found to be compatible. The distributions of the invariant mass m(µ N jet) for same-sign and opposite-sign muon data are shown in Fig. 4 with the fits for the 15 GeV/c 2 neutrino mass hypothesis superimposed. No significant excess of signal is observed in the data for this or any other mass hypothesis. Upper limits at 95% confidence level on B(N → µ jet) |V µN | 2 are set for each heavy-neutrino mass hypothesis using the CLs method [56] with a one-sided profile likelihood ratio [57] as test statistic. The upper limits as a function of heavy-neutrino mass are shown in Fig. 5. For the same-sign muons sample and neutrino mass in excess of 20 GeV/c 2 , the measured limit is between 2 and 3.8 standard deviations above the expected limit. The worse limit obtained with respect to the expectation can be attributed to the four data candidates with m(µ jet) between 20 and 40 GeV/c 2 . The value of the muon identification BDTs for three of the candidates are very close to the requirements, defined a priori with a blinded procedure, indicating that they are background-like and probably a QCD fluctuation. Each candidate has also a relatively large value for the missing transverse momentum in the event, which is not characteristic for the signal. Consequently, the excess at high mass is likely the result of an imperfectly modelled component of the background. For the opposite-sign muons samples, the expected limit is a factor 5 to 10 worse due to the irreducible background from Drell-Yan processes, in agreement with expectations.
To set upper limits on the coupling, the results of Fig. 5 are scaled by B(N → µ jet) = 0.51, computed as described in Section 1 assuming |V eN | 2 = |V τ N | 2 = 0. For the 5 GeV/c 2 heavy-neutrino mass hypothesis, at the limit set, the heavy neutrino is expected to be long-lived with a lifetime of 3.8 ps and 1.1 ps for same-and opposite-sign  muons in the final states, respectively. Since this search targets prompt heavy neutrinos, the acceptance is corrected accordingly. The constraints on the coupling as a function of mass for the opposite-and same-sign muons final state, with and without the acceptance correction factor applied, are illustrated in Fig. 6.

Conclusion
A search for a prompt heavy neutrino in the decay N → µ jet is performed using data from proton-proton collisions recorded by the LHCb experiment, corresponding to a total integrated luminosity of 3 fb −1 . No evidence for heavy neutrinos is observed and limits of the order of 10 −4 and 10 −3 are set as a function of heavy-neutrino mass for leptonnumber-conserving and lepton-number-violating decays, respectively. These represent the first limits on the coupling to a heavy neutrino in the mass range 5-50 GeV/c 2 at LHCb. For the first time the signature of two muons and a low mass jet has been probed for heavy neutrinos with mass lower than 20 GeV/c 2 . Furthermore, this is the first limit on lepton-number-conserving decays of a prompt heavy neutrino in the mass range of interest. The observed limits on lepton-number-violating decays are not yet competitive with the existing limits [4,18,19]. With an integrated luminosity of 50 fb −1 , a better sensitivity than the current most stringent limit could be reached for the same-sign muons channel. While this analysis targets prompt heavy-neutrino decays, better sensitivity for low heavy-neutrino masses can be achieved by including long-lived signatures.