Heavy neutrinos at future linear e+e− colliders

Neutrinos are among the most mysterious particles in nature. Their mass hierarchy and oscillations, as well as their antiparticle properties, are being intensively studied in experiments around the world. Moreover, in many models of physics beyond the Standard Model, the baryon asymmetry or the dark matter density in the Universe are explained by introducing new species of neutrinos. Among others, heavy neutrinos of Dirac or Majorana nature were proposed to solve open questions in High Energy Physics. Such neutrinos with masses above the electroweak (EW) scale could be produced at future linear e+e− colliders, like the Compact LInear Collider (CLIC) or the International Linear Collider (ILC). We studied the possibility of observing decays of heavy Dirac and Majorana neutrinos in the qqℓ final state with ILC running at 500 GeV and 1 TeV, and CLIC at 3 TeV. The analysis is based on the Whizard event generation and fast simulation of detector response with Delphes. Neutrinos with masses from 200 GeV to 3.2 TeV were considered. We estimated the limits on the production cross sections, interpreted them in terms of the neutrino-lepton coupling parameter VℓN2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {V}_{\ell N}^2 $$\end{document} (effectively the neutrino mixing angle) and compared them with current limits coming from the LHC running at 13 TeV, as well as the expected limits from future hadron colliders. The limits for the future lepton colliders, extending down to the coupling values of 10−7− 10−6, are stricter than any other limit estimates published so far.


Introduction
In several models of New Physics, some open problems of the Standard Model (SM), such as the baryon asymmetry in the universe, the flavour puzzle, or the nature of the dark matter (DM), are solved by introducing new species of neutrinos of either Dirac or Majorana nature (see e.g.[1][2][3]).A sector of sterile neutrinos connected to the SM by mixing with the SM neutrinos could exhibit additional CP violation needed to explain the baryon asymmetry in the universe.The lightest sterile neutrino could be stable or so long-lived that it constitutes a considerable amount of the major part of DM.The neutrino sector also plays a prominent role in models with both lepton-flavour violation and lepton-flavour non-universality, which could explain the recent LHC flavour anomalies [4].There are also proposed connections of the anomaly in the magnetic moment of the muon (g − 2) µ to the neutrino sector [5,6].Different mechanisms can be considered for production of such Dirac or Majorana neutrinos with masses exceeding several GeV at existing or future high-energy colliders, e.g.via exchange of SM gauge bosons, additional Z or W bosons, or in decays of heavy new particles like leptoquarks.For scenarios when the production mechanism is via the weak force, lepton colliders seem to be the most suitable devices for the heavy neutrino searches.There are two distinct scenarios."Light heavy neutrinos" with masses below the Z mass can occur in decays of the Z and W boson, and the large luminosity of future Z and electroweak factories (for recent studies cf.[7,8]) would give the best search limits, together with the high-luminosity phase of the LHC (HL-LHC).There is a small intermediate phase where the neutrino would be heavier than W and Z, but lighter than the H(125) Higgs boson.Then, it could occur in (invisible) Higgs decays, but it will be hard to distinguish them from e.g.Higgs portal models.As soon as the neutrino masses are above the electroweak scale, the heavy neutrinos can be produced at future linear e + e − colliders, like the Compact Linear Collider (CLIC) [9] or the International Linear Collider (ILC) [10].The signatures observable at lepton colliders have already been discussed in the literature (see e.g.[11][12][13][14][15][16][17][18][19]), but detailed, quantitative studies taking into account all relevant experimental effects have been missing so far.
Many different heavy neutrino production scenarios have been studied at the LHC.For high masses of the new neutral lepton, above the EW boson masses, the highest sensitivity is expected for the heavy Majorana neutrino searches in the tri-lepton or same-sign di-lepton channels.Limits on the coupling parameter V 2 N extend down to about 10 −5 for neutrino masses between 10 and 50 GeV [20,21], but are significantly weaker for masses above the Z boson mass scale.Limits on the new neutral lepton couplings for masses up to 50 GeV can also be extracted from the analysis of W boson decays [22].Stronger limits, of the order of 10 −6 , were obtained from the search for long-lived particle decays (displaced vertex signature) [21,23], which are however limited to low neutrino masses (below 10-15 GeV).Prospects for heavy Majorana neutrino searches were considered for future hadron colliders [24], as well as electron-proton colliders [25].
In this work, the possibility of observing the production and decays of heavy Dirac and Majorana neutrinos into the qq final state (corresponding to two measured jets and one lepton) at the ILC running at 500 GeV and 1 TeV, and the CLIC at 3 TeV is studied.The analysis is based on Whizard [26,27] event generation and fast simulation of detector response with Delphes [28].Dirac and Majorana neutrinos with masses from 200 GeV to 3.2 TeV are considered.We estimate limits on the production cross section and on the neutrino-lepton coupling using machine learning methods and compare them with current limits coming from the LHC running at 13 TeV, as well as the expected future limits from hadron colliders.Beam-related effects and systematic uncertainties are included in the procedure.The expected limits obtained in our study are stronger than any other estimates published so far and exceed those for pp machines by several orders of magnitude.
The paper is structured as follows: in Section 2, our model setup and simulation framework are described; in Section 3, we present our analysis procedure.Results are discussed in Section 4 and the most important features of the work and prospects of the analysis are summarised in Section 5.

Model setup
There is a vast theory space of models of sterile neutrinos and extended neutrino sectors, regarding which pending problem of the SM is specifically addressed by them: they allow to introduce new sources of CP violation needed for leptogenesis or baryogenesis [29][30][31][32][33], they introduce candidates for (cold) dark matter [34][35][36] and they might play a role in the flavor puzzle [37].Depending on whether they are embedded in extended gauge sectors, like e.g. in left-right symmetric models or Grand Unified Theories (GUTs), there will be additional gauge bosons above the electroweak scale in the multi TeV or not.For this study on the sensitivity reach of future high-energy lepton colliders, we stay mostly modelindependent and assume that -although there are up to three different heavy neutrino The only interaction of the new neutrinos with the SM is through mixing effects, which come from a non-diagonal mass matrix between the electroweak doublet neutrinos and sterile neutrinos.Hence, in this work, we focus on the Phenomenological Type I Seesaw mechanism [38,39], implemented within the HeavyN model [40] with Majorana [41,42] and Dirac [24] neutrinos, an effective extension of the Standard Model introducing three flavours of right-handed neutrinos (denoted as N 1 , N 2 and N 3 ) which are singlets under the SM gauge groups.The Lagrangian of the model is given by: where L N is a sum of kinetic and mass terms for heavy neutrinos (note that we use 4-spinor notation in all cases, which combines terms with spinors of dotted and undotted indices): with an overall factor ξ ν = 1 for the Dirac neutrino and ξ ν = 1 2 for the Majorana neutrino scenarios.L W N corresponds to neutrino interactions with a W boson: L ZN ν to interactions with a Z boson: and L HN ν to interactions with a Higgs boson: The vertices involving N introduced by the model are shown in Figure 1.
The model is described in FeynRules [43,44], the Mathematica package to calculate Feynman rules associated with the Lagrangian of a given model.The output is stored in the UFO format [45], the model format for automatized matrix element generators.The UFO library used in the analysis contains 12 free parameters in addition to the SM parameters: • three masses of the heavy neutrinos: m N 1 , m N 2 and m N 3 , • nine real1 mixing parameters V lk , where l = e, µ, τ and k = N 1, N 2, N 3.
There are also three widths of the heavy neutrinos (Γ N 1 , Γ N 2 and Γ N 3 ) to be set.
For such neutrinos, there are many different signatures expected at future colliders [15].For e + e − collisions, the dominant production channels are s-channel Z production and tchannel W exchange, resulting in the production of a light-heavy neutrino pair: The Z exchange process is dominant at the Z-pole (around the mass of the Z boson), while for centre-of-mass energies above the Z-pole, the W exchange contribution is more important.Analytic calculations show that the cross section for the production of a heavyheavy neutrino pair is much lower and, hence, these processes are not considered in the analysis2 .In the parameter space considered, the heavy neutrino has a microscopic lifetime (cτ 1 nm) so that no displaced vertices can be reconstructed and the products of its decays point back to the primary interaction point.Different final states are possible; however, in this paper, we focus on the qq ν final state, corresponding, at the experimental level, to the jj signature.Example Feynman diagrams for the process are presented in Figure 2. The production process is dominated by the W exchange for which only lefthanded electrons and right-handed positrons contribute and thus, we decided to consider the corresponding beam polarisation settings.Since the signal and the leading SM background channels depend on the polarisation in a similar way, such a choice allows for increasing the expected signal event number, keeping the signal-to-background ratio on the same order.The following collider setups are considered: • ILC500 -ILC running at 500 GeV, with an integrated luminosity of 1.6 ab −1 and beam polarisation of −80% for electrons and +30% for positrons; • ILC1000 -ILC running at 1 TeV, with an integrated luminosity of 3.2 ab −1 and beam polarisation of −80% for electrons and +20% for positrons; • CLIC3000 -CLIC running at 3 TeV, with an integrated luminosity of 4 ab −1 and beam polarisation of −80% for electrons (no polarisation for positrons).
For the dominant production channel, the above runs correspond to about 80% of all data for ILC and 97% for CLIC and the difference is mostly caused by the luminosity fraction assumed to be collected for each polarisation setup at those colliders.

Event generation and benchmark scenarios
The first step was to generate collision events using Whizard [26,27].For the generation of SM backgrounds and Dirac neutrino samples, version 2.8.5 was used, while the simulation of the Majorana neutrino production was the first physics project using the new major version To generate signal events, the Dirac_NLO and Gen3Mass_NLO implementations of the HeavyN model, described within the FeynRules model database, were used.To simplify the analysis, we assumed that only a single heavy neutrino is coupled to the Standard Model particles3 .Therefore, for the simulation, the masses of N 2 and N 3 were set to 10 TeV and their couplings to zero in the model.For the neutrino N 1 that is assumed to have non-vanishing coupling to the SM, and to which we refer from now on just as "heavy neutrino" or N , masses in the range 200-3200 GeV in steps of either 50, 100 or 200 GeV were considered as signal benchmark scenarios.For these scenarios, all the mixing parameters were set equal to √ 0.0003: Widths of the heavy neutrino were calculated using Whizard and are in agreement with the values given in [41].Because of the additional CP-conjugate final states, the widths for the Majorana case are twice as large as for the Dirac case.The width values for the reference scenario are shown in Fig. 3 as a function of the heavy neutrino mass.One can observe that for the assumed coupling, the neutrino can be treated as a very narrow resonance, but the neutrino widths are not so small to produce displaced vertices or even let the neutrinos escape the detector.
As the signal signature, we considered the production of a light-heavy neutrino pair with the heavy neutrino decaying into two quarks (all quarks and antiquarks lighter than t were allowed and their masses were set to zero in Whizard) and one lepton (all flavours allowed, only taus are assumed to be massive), so a decay N → ± jj.For each signal scenario, 300,000 events were generated.The cross section for the process at different collider setups (including beam spectra, beam polarisation and ISR) as a function of the heavy neutrino mass is shown in Figure 4.For masses below the collider energy, the cross section is of the order of 10 fb; then, it decreases fast to 10 For the Majorana case, the width is twice that for the Dirac case (cf.text).also checked that in the wide range of couplings (10 −7 − 1), the cross section can be treated as proportional to V 2 N .For the background samples, the Standard Model implementation (SM ) in Whizard was used, so the processes involving the heavy neutrino are excluded from the background.All the quark, electron and muon masses, as well as the strong coupling constant4 , were set to zero in Whizard to assure consistency with the configuration used for the signal generation.As for the background, we considered processes with at least one lepton in the final state: • e + e − → qq ν, • e + e − → qq ν ν, • e + e − → qqqq ν, • e + e − → qqqq .
Such a choice of background channels was caused by limitations of the detector simulation framework -in Delphes, fake lepton tracks cannot be generated, so at least one lepton in the final state is needed.Events without any leptons would be excluded at the preselection.
Moreover, we included γ-induced background channels.Both beamstrahlung (denoted as B in the following) and photons from collinear initial-state splittings (EPA photon interactions, denoted as E ) were considered in the analysis: • e + γ/γe − → qq (denoted as γe ± → qq ), • γγ → qq ν, where also processes with one beamstrahlung photon and one EPA photon are taken into account.Because of the lack of genuine Circe2 files for the photon spectra of ILC1000, we decided on an approximate solution and scaled the ILC500 spectrum files for usage at a collision energy of 1 TeV, as the shape of the spectra is not expected to change significantly with energy.
One should notice that the expected luminosity for the γ B collisions differs from the e + e − luminosity.The following fractions of the integrated e + e − luminosity are assumed in the analysis: These estimates are based on the detailed simulation of the accelerator performance [48,49].
At the generator level, standard cuts are adopted.We require the invariant mass of the produced quark and lepton pairs to be above 10 GeV and the four-momentum transfer between the outgoing and incoming electrons (or positrons) to be at least 4 GeV.To avoid double-counting, for the EPA events, a maximal photon energy transfer cut of 4 GeV is set.Furthermore, for the samples with beamstrahlung photons, we impose an additional cut on charged leptons to be detected in the central detector (5 • < θ < 175 • , where θ is the lepton polar angle) which helps to remove collinear singularities.
Cross sections for different processes calculated in Whizard are presented in Table 1.

Detector simulation
In the next step, the fast detector simulation framework Delphes [28] was used to simulate the detector response, with cards available for parameterisation of the ILC detector (delphes_card_ILCgen.tcl ) and CLIC detector (delphes_card_CLICdet_Stage3_fcal.tcl ), respectively.As opposed to programs based on full simulation, Delphes provides a general parametrisation of the detector acceptance and response, making the simulation much faster than in the standard approach and allowing for testing many points in the parameter space.
In the ILC detector model, the Durham algorithm was implemented for jet reconstruction, following results of the full simulation studies [50], while for CLIC, the VLC algorithm with the following parameter setup: R = 0.5, β = 1, γ = 1 (see [51] for details) was applied.
Results of the clustering in the exclusive two-jet mode were selected for the presented study based on the expected signal topology.1: Cross section σ and number of expected preselected events N (see Section 3) for different channels at ILC500, ILC1000 and CLIC3000.The cross section for γ B γ B → qqlν at ILC500 is negligible (0.042 fb) because the energy spectrum of the photons is too low for the on-shell W W production.

Analysis procedure
The first step of the analysis was to exclude events resulting in a different topology than the one expected for the signal.Only events consisting of two jets and one lepton (electron or muon) were accepted.Events with any other activity in the detector (additional leptons or reconstructed photons) were rejected.It was also required that the total transverse momentum of final state objects not contributing to the required final state (untagged transverse momentum) had to be smaller than 20 GeV.In particular, this cut rejects events with significant contribution of forward deposits assigned to the beam jets (not included in the final state) by the VLC algorithm.One should notice that events with the qqτ ν final state could also pass the preselection, if the τ decayed into leptons.Numbers of expected events passing the above cuts at the considered future collider options are given in Table 1.
In Figure 5, distributions of the invariant mass of two jets and a lepton are shown for different collider setups.A clear peak corresponding to the heavy neutrino mass is visible in each plot.The left shoulders of those peaks can be explained by the contribution of leptonic τ decays, when two additional escaping neutrinos reduce the invariant mass of the detectable final state.The tails on the right-hand side are caused by detector effects, for example, worse track momentum resolution for leptons going at small angles.It is also important to notice that the background levels for the muon channel are significantly smaller.An extra cut on the invariant mass could help with the background-signal separation at the preselection level, but we do not apply it, as we want to consider broad spectra of heavy neutrino mass values.Nevertheless, it was checked that the cut does not affect the final results obtained with the Boosted Decision Tree algorithm.
In the next step, the Boosted Decision Tree (BDT) method implemented in the TMVA package [52] was used to discriminate between signal and background events.The following 8 variables were considered to train the BDT algorithm: • m qq -invariant mass of the dijet-lepton system, • α -angle between the dijet-system and the lepton, • α qq -angle between the two jets, • E -lepton energy, • E qq -energy of the dijet-lepton system, • p T -lepton transverse momentum, • p T qq -dijet transverse momentum, • p T qq -transverse momentum of the dijet-lepton system.
Other variables were also investigated, but it was found that they did not improve the BDT performance.
The BDT algorithm was trained separately for events with electrons and muons in the final state.The main reason for this approach was the fact that there are more background channels for electrons in the final state and the results for this case were expected to be less stringent.The BDT response for an example reference scenario (Dirac neutrino, m N = 300 GeV) with muons in the final state at ILC500 is shown in Figure 6.In Figure 7, the variable distributions for the same scenario are presented.
In the last step, the CL s method, implemented within the RooStats package [53], was used to extract the cross section limits from the expected BDT response distributions.This approach allows for combining different measurement channels (electrons and muons in this case) and adding systematic uncertainties.We considered only systematics related to the normalisation of the samples.The normalisation uncertainty of the e + e − data sets was set to 1%, and for the γ BS e ± and γ BS γ BS backgrounds additional uncertainties of 2% and 5%, respectively, were applied.These values can be treated as conservative (see e.g.[54]), but it was verified that even without a normalisation constraint (i.e., setting the normalisation uncertainty to 100%), the extracted limits are hardly changed.
We also verified the effect of the jet energy scale uncertainty for a few example points in the parameter space.Jet energy-momentum 4-vectors were scaled up and down by 1%.Since it turned out that there is no impact on the final results, we refrained from studying the effect.Other kinds of uncertainties are also not expected to affect the final conclusions significantly and thus, were not included in the analysis procedure.

Results
After having detailed the analysis methods, we present in this section our results for the sensitivity of ILC and CLIC to heavy neutrinos.In Figure 8, the limits on the cross section   for the considered process are presented, separately for electron and muon channel studies.Better limits for most of the considered scenarios are obtained for muons.Only for the highest neutrino masses at CLIC3000, the limits resulting from the electron channel are slightly stronger.
Extracted cross section limits result from the expected background level and the signal selection efficiency after the optimised event selection.They only reflect the experimental sensitivity and do not depend on the signal cross section predictions.This is why the cross section limits do not get significantly weaker for neutrino masses above the collision energy.Processes mediated by off-shell neutrino exchange are also included in our analysis and signal-background discrimination is only slightly weaker in this case.However, the cross sections for such processes are much smaller than for the on-shell production (refer to Fig. 4), so the corresponding limits on the neutrino coupling V 2 N are much weaker.Such limits are presented in Figure 9, where combined results for Dirac and Majorana neutrino hypotheses are compared.Limits for the two neutrino types are very similar in a wide range of neutrino masses.Below the energy threshold, the differences could be interpreted as statistical fluctuations.However, above the threshold, a separation between the lines is clearly visible.The reason for such a behaviour is the fact that for large neutrino masses, off-shell production above the collider energy is more sensitive to the neutrino width.Since the width of the heavy Dirac neutrino is larger by a factor of 2, so is the production cross section (see Figure 4), and more events are expected to be observed for the same coupling value, resulting in stronger limits.Nevertheless, it has to be noted that the kinematic distributions for the Dirac and Majorana cases are not the same.In Figure 10, the distribution of the lepton emission angle in the N rest frame at the generator level is shown.The flat distribution for the Majorana neutrino corresponds to the isotropic emission (stemming from an average over the two charge-conjugated decay channels), while for the Dirac case, leptons are emitted mostly in the forward direction.In Figure 11, distributions of the lepton energy, E , and reconstructed neutrino mass, m qq , at the detector level, are shown.The lepton energy distribution is significantly affected by the different angular distributions, reflecting the neutrino nature.This was first noticed in [55].However, the distributions of the invariant mass of the qq state are almost identical for both cases.As it is the most efficient criterion for the BDT separation, it explains the very similar results for Dirac and Majorana neutrinos.
The expected limits on the mixing parameter V 2 N compared to current limits and estimates for future hadron machines are presented in Figure 12.The limits for the LHC at 13 TeV come from the CMS Collaboration (Fig. 2 in [20]) and were obtained for neutrinos of Majorana nature, while the limits for future high-energy hadron colliders were taken from the simulation, Fig. 25b, in [24], where Dirac neutrinos were considered.However, when comparing the results, one should note that in the analyses different assumptions regarding the coupling structure have been made: in hadron collider studies, only two non-zero flavour mixings were taken into account, V 2 eN = V 2 µN = V 2 τ N ≡ 0, while all the couplings are assumed to have the same non-zero value, V 2 eN = V 2 µN = V 2 τ N , in our case.Nevertheless, it was verified that our analysis would give even stronger limits if only two non-zero couplings are considered.It is due to the fact that taus can decay into quarks  N for different collider setups (solid lines: ILC500green, ILC1000 -violet, CLIC3000 -dark red).Dashed lines indicate limits from current and future hadron colliders based on [20,24].See the text for details. and then, such events (without electrons or muons in the final state) are excluded from the analysis.On the other hand, as taus can decay leptonically, some of the tau events are included in the analysis, and thus, rerunning of the analysis is needed to compare the results quantitatively with and without employing taus.

Conclusions
Many theories suggest that, beyond the Standard Model, new particles exist.In some models, these particles are neutral leptons with masses above the electroweak scale which could potentially solve observed cosmological problems, such as the baryon asymmetry or the existence of dark matter.One of the ways to search for such heavy neutrinos could be to use future linear lepton colliders.Nowadays, two concepts of linear lepton colliders are considered: the International Linear Collider (ILC) and the Compact Linear Collider (CLIC).For heavy, weak-scale neutrinos, there are a plethora of different models, depending on whether they address primarily the CP problem of baryogenesis, the dark matter paradigm, or whether they are embedded in theories with extended gauge symmetries like e.g.Grand Unified Theories (GUTs).In this paper, we tried to remain relatively agnostic towards specific models and considered only a single kinematically accessible heavy neutrino species, however, allowing for flavour mixing with all three generations.
Neutrinos of both Dirac and Majorana nature and masses in the range of 200 to 3200 GeV were considered.We included all relevant e + e − → X SM background processes, as well as those induced by collinear photon splitting from EPA and induced by

Figure 1 :
Figure 1: Extra vertices in the HeavyN model: interactions of heavy neutrinos with W , Z or Higgs bosons

Figure 2 :
Figure 2: Feynman diagrams for heavy neutrino production at e + e − colliders with the qq ν signature

Figure 3 :
Figure 3: Calculated widths of the heavy Dirac (blue solid line) and Majorana (pink dashed line) neutrinos for the reference scenario (V 2 N = 0.0003).For the Majorana case, the width is twice that for the Dirac case (cf.text).

Figure 4 :
Figure 4: Reference scenario cross section for heavy Dirac neutrino (left plot) and Majorana neutrino (right plot) production, respectively, resulting in the qq ν final state at different collider setups: ILC500 (green solid line), ILC1000 (red dashed line) and CLIC3000 (blue dotted line).Polarisation (left-right for ILC, left-unpolarised for CLIC), ISR photons and beam spectra are included.The reference scenario assumes V 2 N = 0.0003.

Figure 5 :Figure 6 :
Figure5: qq mass distribution for different collider setups for electrons (left) and muons (right) in the final state, respectively.Black solid lines stand for the e + e − background, red dashed lines for the γ-induced background and thick green lines for the different signal scenarios -results for Dirac neutrinos with masses of 300 GeV, 500 GeV and 700 GeV are presented for ILC500, ILC1000 and CLIC3000, respectively.

Figure 7 :
Figure 7: Distributions of the variables used in the BDT procedure for the reference scenario (Dirac neutrino, m N = 300 GeV) with muons at ILC500.The red line stands for the background, the green line for the signal.

Figure 8 :
Figure 8: 95% C.L. limits on the cross section of heavy Dirac neutrino production and decay (the qq ν final state) as a function of the neutrino mass for different collider setups.Dots signify the analysis with an electron in the final state and stars that with a muon, respectively.

Figure 9 :
Figure 9: Comparison between results for Majorana (dashed line) and Dirac (solid line) neutrinos for different collider scenarios: green for ILC500, blue for ILC1000 and red for CLIC3000 respectively.

Figure 10 :Figure 11 :
Figure 10: Distribution of the cosine of the lepton emission angle in the N rest-frame for the Majorana (pink dashed line) and the Dirac (green solid line) neutrinos with a mass of 500 GeV at CLIC3000 (generator level)

Figure 12 :
Figure 12: Limits on the coupling V 2N for different collider setups (solid lines: ILC500green, ILC1000 -violet, CLIC3000 -dark red).Dashed lines indicate limits from current and future hadron colliders based on[20,24].See the text for details.