Collider Searches for Heavy Neutral Leptons: beyond simplified scenarios

With very few exceptions, the large amount of available experimental bounds on heavy neutral leptons - HNL - have been derived relying on the assumption of the existence of a single (usually Majorana) sterile fermion state that mixes with only one lepton flavour. However, most of the extensions of the Standard Model involving sterile fermions predict the existence of several HNLs, with complex mixing patterns to all flavours. Consequently, most of the experimental bounds for HNLs need to be recast before being applied to a generic scenario. In this work, we focus on LHC searches of heavy neutral leptons and discuss how to reinterpret the available bounds when it comes to consider mixings to all active flavours, not only in the case with a single HNL, but also in the case when more heavy neutral leptons are involved. In the latter case, we also consider the possibility of interference effects and show how the bounds on the parameter space should be recast.


Introduction
Generating neutrino masses and their mixing as observed in neutrino oscillation phenomena requires to go for beyond the Standard Model of Particles (BSM). Many options are presently explored as extensions of the Higgs and/or gauge sectors, most of the time with new fields within the particle content. In particular heavy neutral fermions, such as right-handed neutrinos ν R , are often present as building blocks of several neutrino mass generation mechanisms. For instance, at least two ν R are required to accommodate light neutrino masses via the type-I seesaw mechanism [1][2][3][4][5]. Moreover, in several variants of the type-I seesaw realized at low scale, other sterile (from SM gauge interactions) fermions ν S are considered, as in the case for the Inverse [6][7][8] and Linear [9,10] seesaw mechanisms; these variants allow to have large neutrino Yukawa couplings with a comparatively low seesaw scale, potentially within reach at colliders. From now on, we will refer to these states (ν R,S ) as Heavy Neutral Leptons (HNL).
Due to the presence of HNLs, the charged and neutral currents are modified, with the leptonic mixing matrix encoding now not only the PMNS mixing matrix [11,12], but also the active-HNL mixings U αN , α = e, µ, τ . With these modifications and depending on the mass scale of these new neutral leptons, one expects an impact on numerous observables and thus to abundant constraints on the plane (M N , |U αN | 2 ) (see Refs. [13][14][15] and references therein).
In this work we will be interested in HNL searches at high-energy colliders, mainly at the LHC due to its current extensive program dedicated to HNL searches, see for instance [16], and the numerous dedicated works and analyses . With very few exceptions, the large amount of available HNL bounds have been derived relying on the assumption of the existence of a single (usually Majorana) HNL that mixes with only one lepton flavour. However, most of the BSM scenarios involving new neutral leptons address the lepton mixing as a whole, as it impacts flavour physics studies and lepton properties -Dirac or Majorana nature for neutrinos -with different lepton number conserving/violating processes. The mixing pattern in these scenarios is expected to be quite complex, so applying the bounds from negative searches on the HNL space of parameters derived in the context of simplified hypotheses (with only one HNL which mixes to only one active flavour) does not seem adequate. Indeed, as we will see, using these limits directly will in general overconstrain the parameter space. Consequently, most of the experimental bounds for HNL need to be recast before being applied to a generic BSM scenario.
The motivation for reinterpreting LHC bounds in general is a well-established topic, see for instance Ref. [59] and references therein. In the context of HNL searches, the reinterpretation of the obtained bounds on the HNL mixings to active flavours has been addressed in previous works, as in for instance Ref. [50] discussing searches for Heavy Neutral Leptons with Displaced Vertices, and more recently in Ref. [58] focusing on HNLs decaying promptly to a tri-lepton final state. In the latter, the single-flavour mixing results obtained by ATLAS have been recast to a low-scale seesaw model with a pseudo-degenerate pair of HNLs, the most minimal and simple extension in order to accommodate neutrino oscillation data (the lightest neutrino being massless). Due to the simplicity of this model, the active neutrino masses and mixings determine the flavour pattern of the HNLs [60,61], and it is possible to define benchmark points beyond the single-flavour scenario [62]. While being a very interesting scenario, as it is probing the parameter space connected to light neutrino masses, this approach has the drawback of being model dependent. For example, considering other sources for light neutrino masses 1 , such as additional HNLs not necessarily within the LHC range, could spoil the correlation between light and heavy sectors that motivated the definition of these scenarios.
For this reason, in this work we will follow a different approach. We will instead work with physical HNL states with independent mixings and masses, with the motivation of covering every scenario that could be realized at generic BSM models. We will also discuss how to go beyond the simplest single-flavour mixing scenario, however we will not attempt to explain light neutrino masses and mixings. This idea is actually the most straightforward extension to what is usually assumed at LHC searches. In doing so, we will discuss what would be the most relevant quantities to be bounded experimentally in order to easily reinterpret the results.
Furthermore, we will also extend the study to the case where more than one HNL is present and, in particular, taking into account possible interference effect when at least two heavy neutral leptons are in the same mass regime (nearly degenerate or possibly forming a pseudo-Dirac neutrino pair). This possible interference effect can lead to different bounds on the active-sterile mixings, see for instance [63][64][65].
This work is organised as follows: after having thoroughly discussed the status of high-energy collider searches of HNL in Section 2, a first insight of going beyond the single active-sterile mixing approximation is examined in Section 3. In this section, we also discuss how one would reinterpret experimental data from branching ratios to limits on the model parameter space, i.e. HNL masses and mixings. In Section 4, we consider the generic case of having more than one HNL and discuss possible changes in the interpretation of the bounds on HNL, in particular when the interference effects can be relevant. Final comments are collected in Section 5. Further details for deriving the relevant amplitudes are collected in the Appendices.

Status of HNL searches at high-energy colliders
Heavy Neutral Leptons can be searched for in a wide variety of processes and experiments, the HNL mass being the key ingredient to decide which is the optimal one. HNLs lighter than the GeV scale can lead to signatures in nuclear β decays or in leptonic or semileptonic meson and tau decays, while heavy HNL above the TeV are better explored indirectly by electroweak (EW) precision observables or rare flavour processes. For a detailed review of all these signals and experimental status, see for instance Refs. [13][14][15].
Here, we are interested in the intermediate regime, with HNL masses M N ranging from few to hundreds of GeVs. Such HNLs could be directly produced at high-energy colliders and their lifetimes are usually short enough to decay within the detectors, so we could discover them looking for their decay products. As these are weak processes, the Standard Model boson masses obviously define the relevant scale to be compared to. Along this work, we will refer as light 2 HNL to those with masses lighter than the W boson, and as heavy HNL to the ones with M N > M W . Extensive reviews about HNL searches at colliders can be found for instance in Refs. [66][67][68]. Here we just summarize and update the list of experimental analyses, introducing at the same time the most relevant aspects that we will use in our discussion in the next sections.
As in any collider search looking for heavy unstable particles, we need to consider both the production and the decay channels of the HNLs. At a hadronic collider such as the LHC, the main production channel comes from Drell-Yan W and Z bosons, where the gauge bosons could be on-or off-shell, depending on whether the HNL is lighter or heavier than the W or Z bosons. Additional production channels could also arise from the Higgs (H) boson decays, which could be motivated in several models providing large neutrino Yukawa couplings. Unfortunately, Higgs bosons are produced less abundantly than weak bosons, so they are usually neglected. Moreover the W channel has the additional prompt charged lepton that can help triggering the process and reducing backgrounds, and thus experimental searches focused mostly on this channel. Nevertheless, it is worth mentioning that for very heavy masses, at around the TeV scale, vector boson fusion channels such as W γ or W W become important and could even dominate the production of HNLs [25,27,57]. Indeed, the latest CMS analysis [69] already included the W γ channel in order to enhance their sensitivity to high HNL masses. After being produced, a HNL of several GeVs, but still lighter than M W , will decay dominantly via off-shell W or Z bosons to a 3-body final state On the other hand, if M N is above the EW scale, the dominant decays will be to on-shell W, Z and H bosons, i.e. N → ± W ∓ , νZ, νH. These 2-body decays will be followed by the decay of the heavy bosons, leading at the end to the same final states as before. Nevertheless, it is important to keep in mind that the kinematics in these two mass regimes will be different. Combining both production and decay channels, we get a full process such as the example shown in Fig. 1. Depending on the relative size of M N and M W , either the first or the second W boson will be on-/off-shell, distinguishing the light and heavy HNL regimes. A full catalogue of HNL signatures, combining the different production and decay processes, can be found for instance in Ref. [70]. Here we focus only on those that have been already searched for at the LHC, which we collect in Table 1. We notice that similar searches considering the existence of right-handed currents have also been performed [71][72][73][74][75][76][77][78], and could in principle be recast to our setup. Nevertheless, one would naively expect lower sensitivities, as they are optimized for heavy right-handed gauge bosons.
Most of the LHC searches focused on the smoking gun signature for Majorana neutrinos, the same sign (SS) dilepton final state: W W q q ℓ α ℓ β q q N Figure 1: Drell-Yann HNL production leading to a dilepton signature. The thunder-shaped arrow indicates that the HNL could be of Dirac or of Majorana nature, short-or long-lived.
Here, the lepton pair is accompanied by at least two jets (see Fig. 1), unless M N is much lighter or much heavier than M W , which leads to boosted objects and collimated jets that are reconstructed as a single one. Being a LNV process, the SS dilepton does not suffer from severe SM backgrounds. Unfortunately, current collider searches are sensitive only to relatively large mixings between the HNL and the active neutrinos, too large to explain the masses of the light neutrinos unless a symmetry protected scenario is invoked. More specifically, this symmetry is an approximated conservation of lepton number [79], which also suppressed the expected LNV signal from HNLs (see however Refs. [80][81][82]).
From this point of view, searching for opposite sign (OS) dileptons, as done by LHCb [83], seems more relevant to explore theoretically motivated scenarios. The drawback is the large amount of background from Z → + − decays, which reduces the sensitivity. A possible alternative would be focusing on LFV channels to reduce backgrounds [32,52].
Yet another alternative considers the fully leptonic process This channel has a trilepton signature, rather clean in a hadronic collider, nevertheless it also has a source of MET, which might spoil the complete reconstruction of M N . The trilepton channel offers the possibility to search for both LNV and LNC signals, however most of the experimental analyses still focus only on the LNV channels as to reduce backgrounds, again from Z → + − . For example, ATLAS searched [84] for e ± e ± µ ∓ and µ ± µ ± e ∓ channels 3 , but not for e ± e ∓ µ ± and µ ± µ ∓ e ± . CMS did something similar for light HNLs, although they also included channels with OS but same flavour lepton pairs in the heavy HNL regime, removing only those events with lepton pairs compatible with a decay of a Z boson [69]. As stated before for the SS dilepton channel, searching for HNLs without assuming their Majorana nature will be helpful to probe scenarios compatible with neutrino oscillation data, and thus with potentially suppressed LNV signals. Finally, it is important to stress that improving the experimental sensitivities to smaller values of mixings implies exploring HNL with longer lifetimes, which can travel macroscopic  distances before decaying. Such long-lived HNL would avoid the searches mentioned so far, as they all assumed prompt decaying HNLs, and therefore we need a dedicated search for this kind of topologies. Recently, both ATLAS [85] and CMS [86] have searched for displaced vertex signatures for light HNLs, setting the strongest constraints for GeV masses up to 20 GeV. This can be seen in Fig. 2, where we summarize all the relevant LHC constraints explained in this section.

Searches at lepton colliders
Despite the fact that our current most powerful high-energy collider is a hadronic one, it is important to stress that lepton colliders are extremely relevant for HNL searches. Not only due to the impressive sensitivities expected at future leptonic colliders such as the FCCee [93], but also because HNL searches at LEP still provide the most relevant limits for some M N hypotheses. The great advantage of a leptonic collider is its clean environment, in contrast with the hadronic ones. In the case of LEP, they combined this cleanliness with the huge amount of Z bosons they collected to search for HNLs produced in Z → νN . Moreover, they considered both visibles and semi-invisible HNL decays, such as monojet final states [94]: with the qq pair clustered as a single jet due to the large HNL boost (efficient for M N 30 GeV). For heavier masses, M N ∈ (30, 80) GeV, the signature was composed by two jets with or without a charged lepton. Such a search would be very challenging at a hadronic collider, however it has the advantage of being sensitive to all flavours, including the mixing to the τ lepton, not explored so far by LHC searches. The DELPHI resutls [94], derived for both long-lived and prompt light HNLs, were not improved (for mixings to e and µ flavours) by LHC until very recently, and still dominate for some mass ranges (cf. Fig. 2). Additionally, the L3 collaboration explored the heavy HNL regime by considering their production via the t-channel W diagram [95]. This process dominates the heavy HNL production  Table 1. In the upper (lower) panel a single mixing scenario to electrons (muons) is assumed. Shadowed area cover the area excluded by direct searches at LEP. Notice that below 2 GeV and above (approx.) 100 GeV, bounds from meson decays and from non-unitarity of the lepton mixing [92] dominate respectively over current LHC bounds, although we do not show them explicitly for easier reading of the collider results.
at a e + e − collider running above the Z pole, however it is sensitive only to mixings to electrons. The results by L3 still provide the strongest limits for masses between 100 and 200 GeV. Despite the great effort in the search for HNLs by both LEP and LHC, it is important to analyze their implications for realistic models introducing and motivating the existence of HNLs. A common feature of all these searches is the assumption of a simplified scenario, most of the time consisting on a single HNL mixing to a single lepton flavour, which is not the standard hypothesis one would use from the theory side. To our knowledge, the only exceptions to these simplifications are provided by the recent ATLAS search for long-lived HNLs that also considered a minimal but realistic 2HNL scenario [85], and CMS searches for SS eµ final states [69,89], although still neglecting the mixing to taus. In the two following sections, we discuss the importance of going beyond these simplified scenarios to explore more realistic scenarios.

Beyond the single mixing assumption
As explained in the previous section, most of the LHC analyses are done assuming the existence of just one HNL that mixes to a single flavour, which we referred to as the single mixing scenario. In this section, we consider deviations from this simplified hypothesis and discuss their implications for reinterpreting the LHC bounds summarized in Fig. 2. In particular, we focus only on prompt searches, while the implications for long-lived HNL were discussed in, for instance, Ref. [50].
For simplicity, we still consider the presence of a single HNL (we will discuss deviations from this hypothesis in the next section), however we open the room for generic mixing patterns. Moreover, we will follow a bottom-up approach where the SM is extended by ad-hoc masses for the 3 active neutrinos, as required by oscillation phenomena, and by the presence of the additional HNL N . This framework is useful to study, in particular, the collider phenomenology of HNLs without assuming any specific underlying model or mechanism of light neutrino mass and leptonic mixing generation 4 .
In such a framework, the lepton mixing matrix is thus enlarged to a 4 × 4 unitary matrix so the would-be-PMNS matrix U νν is no longer a unitary matrix, a feature which is indeed used to constrain these models [92], and the fourth column contains the HNL mixings to each flavour: For our discussion, it is interesting to parametrize this column as where U 2 represents the total (squared) mixing of the HNL and the ε α its flavour strengths, with |ε e | 2 + |ε µ | 2 + |ε τ | 2 = 1. Notice that this framework is precisely the one considered by most LHC analyses, the only difference being that they simplify it by setting the a priori non-relevant mixings to zero. Here, we are interested in knowing how these bounds need to be modified in a generic mixing pattern scenario. The reason why we expect the bounds to be modified is twofold. The first reason is due to the importance of the HNL decay width, which depends on every mixing U αN , and which plays a major role in the resonant searches (on-shell produced HNL) we are interested in. This means that the final cross sections will depend on all of the mixings, even on those flavours that are not explicitly present in the charged leptons involved in the process. The second reason is that for some channels, considering generic mixings could open new contributing diagrams, which could modify the distributions and thus the efficiencies of the searches as discussed thoroughly in Ref. [58]. Complete expressions for the computation of the total decay width Γ N can be found for instance in Ref. [96], however for our purposes we parameterize it as where Γ α N stands for the sum of partial decay widths depending on the mixing U αN , after factorizing the |U αN | 2 dependence itself. Thus, Γ α N are independent of the mixings (at leading order) and depend only on the HNL mass. Moreover, when the HNL is heavy enough so that we can neglect charged lepton masses, we get Γ e N Γ µ N Γ τ N and thus With this discussion in mind, we can now study how the different processes displayed in Table 1 depend on the HNL mixings. Let us start focusing on the dilepton channels, the most straighforward case as we only need to track the effect of the HNL total decay width.
In the narrow width approximation, the processes with SS and same flavour dileptons can be factorized in the production of the HNL together with a charged lepton, times its subsequent decay to the same lepton plus jets. The first part depends only on the mixing to the flavour of that lepton, however the second one involves all the mixings due to the HNL decay width. More explicitly, we have or, assuming a heavy enough HNL, Then, we clearly see that those bounds obtained in the single mixing benchmark (ε α = 1) will be relaxed in a general flavour scenario with a fixed U 2 , since in general we will have |ε α | 2 ≤ 1. This is actually the expected behaviour, as switching on other mixings opens for new decay channels, so not every produced HNL will decay to the final state we are searching for. We can repeat the exercise for the different flavour SS dilepton processes. Obviously, the minimal setup in this case requires to have two non-zero mixings, which leads to two diagrams that in principle interfere. Nevertheless, using the narrow width approximation, we can see that both diagrams cannot resonate at the same time, so we can neglect the interference and add both processes incoherently: The case of the trilepton channels can be more involved, mainly because we cannot know the lepton number and the flavour carried by the missing (anti)neutrino. Let us consider first the case of same flavour trileptons. There are two contributing diagrams, one with a neutrino and one with an antineutrino, which we can add incoherently 5 . Then, considering for simplicity a W + Drell-Yan channel, we have process (prompt) Relevant parameters (Majorana HNL) approx.
complete dependence When the trilepton signal involves leptons of two different flavours, we need to consider two subcases: the same-sign same-flavour (SSSF) and opposite-sign same-flavour (OSSF). These channels are trickier because in a generic flavour pattern there are new diagrams not present in the single mixing scenario. For instance, in the case of the SSSF, we have two types of contributions: It is clear that the first process requires of a Majorana HNL, while the second one needs mixings to both flavours. As before, the interference is negligible, so we have Due to the missing (anti)neutrino, both processes are almost identical at the LHC, with the only difference coming again from the different distributions and acceptances of the experimental analysis. This was studied in detail in Ref. [58] for the case of light HNLs. Nevertheless, in order to get a first rough estimate, we can neglect these differences and write The case of OSSF is similar, although now the roles of Majorana and Dirac HNLs are flipped: process (prompt) Relevant parameters (Dirac HNL) approx.
complete dependence Moreover, since we are working in the prompt HNL regime, we can also have pp → + β N, N → ν α . This means that a proper recasting of this kind of signals would require to compute the efficiencies for all these diagrams. If, for the shake of this discussion, we neglect these effects again, we get: Finally, and even if no LHC searches have been performed so far, we can also consider the case with 3 different flavours. Following the same steps, we get We summarize our discussion in Table 2. Here, we assume a Majorana HNL, although a similar table can be easily obtained for Dirac HNL by just switching off the LNV channels we discussed above. This is done in Table 3. Notice how some channels that were designed for Majorana HNLs are also sensitive to Dirac HNLs in the case of having generic mixing patterns, as it was already discussed in Ref. [58].
These tables are to be compared with the minimal mixing scenario where all of the processes scale as |U αN | 2 , or |U αN | 2 |U βN | 2 /(|U αN | 2 + |U βN | 2 ) for α = β. However, we see that in general each process is sensitive to a different combinations of mixing strengths. This means that, in order to generalize the bounds to a generic pattern, it is better to set bounds on the quantity on the last columns of Tables 2 and 3, since then we only need to recompute the new BRs for each mixing hypothesis.
While there are some experimental results also presenting (in the single mixing assumption) bounds in the (M N , |U αN | 2 ×BR) plane, most of the results are given directly in the (M N , |U αN | 2 ) one. Translating the latter to the former is straightforward in most of the cases, since the experimental collaborations usually assume a constant 6 BR for channel under study, although not always specifying the precise value they used. Nevertheless, this recasting to |U αN | 2 ×BR is not possible when the experimental results on |U αN | 2 are presented after combining different channels. This is the case for instance for the latest CMS searches for trileptons [91], where they combined channels like e ± e ± e ∓ and e ± e ± µ ∓ (see Table 1). While in the single mixing scenario both channels depend only on |U eN | 2 , they have a different dependence in the case of a generic mixing scenario (see Table 2), and thus it is not easy to recast the obtained bounds without a dedicated analysis. For this reason, together with the potential efficiency differences discussed above, we will focus the rest of our discussion only on the dilepton channels. As a case of study for the use of Table 2, let us consider the CMS [69] and LHCb [83] searches for the LNV dimuon channel µ ± µ ± . For a given total active-sterile mixing U 2 , we can display the full flavour mixing space in a ternary diagram, as in the left panel of Fig. 3. Then, the single mixing scenario constraints by both CMS and LHCb lie in the top corner. As we move along the ternary, we decrease the flavour strength to muons, so the bounds are relaxed, as shown in the right panel (dashed lines). Here we chose just few benchmark points for the light HNL mass regime, although the same logic applies to the heavy one.
As discussed above, the physical reason for the relaxation of these bounds is due to the new HNL decay channels in the generic mixing scenario. On the other hand, this also implies that the experimental searches for the other channels with different flavours might become relevant. In order to show the interplay between different flavour channels, let us consider again a ternary diagram. We can understand it as a subspace of the parameter space with fixed values of M N and U 2 , which is dissected in flavour space. Then, searches for dimuon channels will cover the area close to the ε µ = 1 corner, becoming weaker as we move further away. Equivalently, dielectron searches will cover the ternary from the ε e = 1 corner, while the e ± µ ± channel will cover the area in between. This is depicted in Fig. 4 for two benchmark points, one in the light  (Table 1), which were derived within a single mixing scenario (or assuming U τ N for the e ± µ ± channel). The white area is still allowed by LHC searches. regime (M N = 30 GeV) and one in the heavy regime (M N = 300 GeV). Figure 4 clearly shows the complementarity of the different dilepton searches, to which we could supplementary add the bounds from trilepton channels if the above mentioned concerns are solved. If the combination of every channel covered all the area of the ternary, we could say that this (M N , U 2 ) point is excluded no matter which flavour mixing pattern we were considering. This is not the case in any of the two examples in the figure, since the bottom left corner is still allowed by LHC. Notice that this is always the case at present, since there are no LHC searches for HNLs mixing to the tau lepton, so the corner of ε τ = 1 will always be allowed. Although we already discussed that a single mixing scenario does not seem very natural for a realistic model including HNLs, closing this gap still motivates the need of performing dedicated searches in the tau sector.
On the other hand, it is important to stress that the benchmark points in Fig. 4 are chosen just for illustrative purposes, since they are already excluded by LEP searches or global fit bounds. Indeed, this is actually the case in a large part of the parameter space for the current LHC bounds discussed in Sec. 2. Nevertheless, it is worth emphasizing that LHC sensitivities are expected to improve during the currently ongoing runs, pushing our knowledge about HNLs beyond present limits. Therefore, combining the different LHC channels as we discussed here will become crucial in order to determine whether a heavy neutral lepton with a given mass and mixing is completely excluded or not.

Beyond the single HNL
Should HNL exist in Nature, there is a priori no bound on their number and in general, BSM models involving HNLs do not introduce just one of them. For example, in the standard type-I seesaw model at least two HNL are needed to explain neutrino oscillation data. In this work, in order to explore deviations from the single HNL hypothesis, we extend the framework considered in the previous section to include more neutral fermions focusing here on the minimal case of SM extension via two HNLs. In this case, the lepton mixing matrix is now a 5×5 unitary matrix with the fourth and fifth columns encoding the mixings of both HNLs to the active leptons For the sake of simplicity, we will assume in this section that the two HNLs mix to just a single flavour. In a sense, it can be seen as extending the interpretation of the bounds of Fig. 2 in horizontal in U νN , to additional columns, while in Section 3 we extended them in vertical, to additional flavours. The most general case with several HNLs and arbitrary mixing patterns can then be inferred as a combination of these two discussions.
In the case where only one of these HNLs is within experimental reach or when there are several HNLs but with well-separated mass regimes, our conclusions derived following the single HNL scenario will apply to each of the HNLs. Nevertheless, if two HNLs happen to be close in mass (as motivated by low-scale seesaw models [61,97,98], resonant leptogenesis [99] or ARS leptogenesis [100]), they could lead to interference effects and modify the results and the bounds obtained in the single HNL hypothesis. Moreover, these modifications might affect both LNV and LNC branching ratios, and thus studying their correlation could shed light on the nature of the HNL (see for instance [64] and references therein).
In this section we discuss how these effects could affect the LHC bounds obtained in the single HNL scenario from searches for LNV and LNC channels, and provide a recipe to combine both results (on LNV and LNC searches) in order to bring forth more robust bounds on HNLs parameter space. The recipe is also applicable to the case where HNLs interfere. We mostly focus on the LHCb [83] results for the prompt dimuon channel, since this is the only available analysis addressing both SS and OS dilepton channels (cf. Table 1). More specifically, this search considers light HNLs which are dominantly produced from on-shell DY W bosons in Fig. 1, that is, In presence of just a single Majorana HNL, although the angular distributions will be different (see Appendix A for more details), the predictions for the total rates of these two processes are of similar size: equal for channels with α = β and twice as large for α = β. The reason for this difference is that the channel with crossed α and β is also contributing in the case of the LNV process. However, it must be added incoherently to the process since the rate is dominated by on-shell HNLs, fixing the momentum of the first lepton with the 2-body decay kinematics, and thus the interference becomes subdominant. This means that for channels with α = β, the total rate for the LNV process is enhanced by a factor of 2 with respect to the LNC one. On the other hand, in channels with α = β there is an additional 1/2 factor from having two identical particles, and thus we obtain the same total rate for both LNV and LNC. When assuming the existence of two HNLs, we have two identical contributions to the total amplitude of the processes, one for each N i . The squared amplitude is then given as sum of the individual contributions of each HNL, plus a potential interference between N 1 and N 2 contributions. For each individual contribution to the amplitudes, we find that they are proportional to U * αN i U * βN i for the SS process and to U * αN i U βN i for the OS one, proving convenient to define with φ αi ∈ [0, 2π], α = e, µ, τ and i = 1, 2. In this way, the interference term resulting from both amplitudes will be proportional to where we have defined δφ ± as follows: with +/− for the SS/OS channel. Details with the complete analytical amplitude and decay rates, which we have additionally checked with Whizard [101], can be found in Appendix A, both for the case with just one HNL and with two HNLs. As already stated, when the mass difference of the two HNLs is too large (compared to their decay width), the contribution of each HNL resonate independently and the interference is negligible. In this case, the total rates for the LNV and LNC rates are related as in the case of the single HNL scenario. Therefore, in the following we will assume the scenario when the interference effects can modify considerably the relative predictions between LNV and LNC rates, which corresponds to the case where both HNLs are close in mass. More specifically, we assume that the individual contribution of each HNL is of similar size, However, we consider that the two HNL mass splitting ∆M N ≡ M N 2 − M N 1 could be different from zero as long as it is small compared to the decay width Γ N . Notice that this kind of scenario appears naturally in low-scale seesaws due to the approximated lepton number conservation. The only difference is that in these models the phases of the HNLs are fixed to be opposite, i.e. δφ + = π and δφ − = 0 (the heavy neutrinos forming a pseudo-Dirac neutrino pair), and that ∆M N is somehow related to light neutrino masses [81,82], while in the following we will let them free, in the spirit of the bottom-up approach described before.
Under these conditions, we can write the total decay rate driven by the two HNLs in the case of W + channel as and equivalently for the W − channel with K − (y, δφ ± ) . Here, we have factorized the total rate in presence of only one HNL, and defined the modulation function K ± y, δφ ± ≡ 1 + cos δφ ± 1 1 + y 2 ∓ sin δφ ± y 1 + y 2 , The function K ± codifies the role of the interference. In the limit of ∆M N Γ N , the two HNL are too separated in mass, coherence is lost and the total contribution is just twice the single HNL contribution for both LNV and LNC. On the other hand, for ∆M N < Γ N , the modulation function can take values from 0 (maximally destructive interference) to 2 (maximally constructive). Thus, we are maximizing the effects of the interference between the two HNLs. Moreover, these effects will be different for LNV and LNC, breaking the equal size prediction in the single HNL scenario.
These modulation function can be used to simply recast the bounds derived by LHCb [83] under the assumption of a single HNL. Noticing that LHCb searched only for the W + channel and given Eq. (15), to recast these bounds to our scenario with two HNLs, we need to rescale the mixing as with δφ − = 0 in this channel with α = β. Following this modulation, we show in the left panel of Fig. 5 how the LHCb bounds might vary, depending on the values of y and the relative phases δφ + , for both LNV and LNC searches, which defines the green and blue bands, respectively. The vertical axis needs to be understood now as the (squared) mixing of each of the HNLs to muons, which in absence of interference (y 1) is just a factor of two stronger with respect to the single HNL scenario. For y 1, however, constructive interference can strengthen the bounds by up to an additional factor of two, while destructive interference could relax it, even avoid it completely in the case of LNV signals. The latter corresponds to the case where δφ + = π, which is precisely when the two HNLs have opposite phases (thus forming a pseudo-Dirac pair), as required by low-scale seesaws with approximated lepton number conservation. Interestingly, the same choice of parameters that maximizes the destructive interference of LNV channels also maximizes the constructive interference for the LNC ones, making the bounds from the latter stronger.
This interplay prompts us to consider both channels at the same time, using their complementarity to set absolute bounds on the mixings that could not be avoided even with ad-hoc values of the parameters that maximize the interference. We show this in the right panel of Fig. 5 for a particular example of M N = 30 GeV and opposite phases, δφ + = π, as required by low-scale seesaws. When the contributions of the two HNLs is maximally coherent (y 1), the LNV searches are avoided at the price of maximizing the LNC bounds. If coherence is lost (y 1), then the stronger LNV bounds always dominate over the LNC ones. We see then that the largest possible mixing could be obtained in between the two cases, when the two tendencies cross over, which we can consider as an absolute bound on the mixing that cannot be avoided even with 2 interfering HNLs.
To summarize, searches for LNC and LNV processes are important and complementary, since they can cover areas of the parameter space even in the case where there is some interference (partially) cancelling any of the two channels. This strongly motivates the need of searching for both LNC and LNV channels in parallel, even if the latter is more challenging experimentally, since combining both of them we could set more robust bounds on generic scenarios including more than one heavy neutral lepton.

CP violation
When a quasi-degenerated pair of HNLs is considered, new CP-violating phases are introduced, which can induce differences in the decays of these particles to leptons over antileptons. If HNLs were discovered in processes such as those in Eq. (28), and provided that enough events were collected, one way of measuring this potential CP asymmetry would be by defining the ratio [65] Using our previous results for the W decays, see Eq. (32), it is straightforward to see that A ± CP takes the simple form where +/− denotes again LNV/LNC processes. As can be seen, this equation does not depend on the HNL masses, but on their mass difference ∆M N (through y = ∆M N /Γ N ). It vanishes for the obvious case of δφ ± = 0, since then there is no difference in both the W decay and its CP equivalent; notice that it is also the case when y is close to zero or too large. This is shown in Fig. 6 where we display the CP asymmetry for different values of δφ ± .
Notice that a similar computation was performed in Ref. [65] with which our results agree but with some slight differences. On the one hand, we find a discrepancy when applying the narrow width approximation for the interference term (see App. A for more details). On the other hand, we performed the computation considering very prompt heavy neutral leptons, so we did not take into account their time evolution and their possible oscillations before decaying. Our approach is appropriate for heavy HNLs, while one should take into account the evolution effect for very light HNLs, thus with longer decay lifetimes, as it was done in Ref. [65]. Finally, we also point out that, in order to measure A CP at a proton-proton collider, we must take into account that the production rates for W + and W − are not the same [65].

Conclusions
In this work, we focused on LHC searches of heavy neutral leptons that decay promptly (shortlived). In most of the searches, which we summarized in Sec. 2, the obtained bounds were derived under the hypothesis of the existence of a single (usually Majorana) HNL that mixes with only one lepton flavour, while most of the BSM scenarios involving new neutral fermions require more than one HNL. Moreover, unless some specific symmetries are present, the mixing pattern in these BSM scenarios is more complex, each HNL mixing in general with all charged leptons, and thus, the bounds derived from negative searches on the HNL parameter space have to be recast before being applied to a generic BSM scenario.
In this study we discussed how to recast the present experimentally obtained bounds on the parameter space, i.e. active-sterile mixings U αN , α = e, µ, τ , versus the mass of the HNL, to the case of generic mixing to all active flavours as well as to the case with several HNLs. The former was covered in Sec. 3, where we inspected the flavour dependencies of each of the channels searched for by the LHC, and stressed the importance of setting bounds not only on the mixings, but on the relevant combination of |U αN | 2 × BR (see Tables 2 and 3). Considering the bounds on this combination, we proposed a method to combine the results in flavour space, using the ternary diagrams in Fig. 4 to conclude whether that area of parameter space is fully excluded, regardless of the assumed mixing pattern.
In the case with several heavy neutral leptons, we focused on the scenario when two HNLs are in the same mass regime, nearly degenerate or possibly forming a pseudo-Dirac neutrino pair, paying special attention to the non-trivial role of the interference between their contributions. To illustrate its importance, we focused on dilepton channels and moreover considered both channels with same and opposite charge of the final leptons, as it was done by LHCb [83]. We showed the complementary of the LNC and LNV searches and the importance of performing both of them in parallel. We stressed that by doing this, we are not only taking into account the two possible nature of a single HNL, Dirac or Majorana, but also covering the case when, for example, two Majorana HNLs exist and however they interfere destructively, suppressing the expected LNV signature.
To summarize, we have discussed the importance of going beyond simplified scenarios such as the single mixing hypothesis. While they are useful for simplifying experimental analyses, they are not directly applicable to BSM models introducing HNLs. Unfortunately, recasting the bounds of each experimental analysis to a given BSM scenario can be a tedious task. In this work, we have proposed an alternative way of presenting the bounds on the parameter space of the HNLs, which under some approximations can be both directly constrained by experimental analyses and also easily recast to a generic BSM scenario. Acknowledgements PE is grateful for the hospitality of the Pôle de Physique Théorique of the IJCLAB (Orsay) during the first stages of the project. This project has received support from the European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie grant agreement No. 860881 (HIDDeν network). The work of PE is supported by the Spanish grants PID2020-113775GB-I00 (AEI / 10.13039/501100011033), CIPROM/2021/054 (Generalitat Valenciana) and by the FPI grant PRE2018-084599. XM also acknowledges partial support from the Spanish Research Agency (Agencia Estatal de Investigación) through the Grant IFT Centro de Excelencia Severo Ochoa No CEX2020-001007-S and Grant PID2019-108892RB-I00 funded by MCIN/AEI/10.13039/501100011033.

A Decay process in presence of 2 HNL
In this appendix we collect the relevant details for the computation of the W boson decay into a HNL, followed by its semileptonic decay. This is the relevant process for the searches performed by LHCb [83] and that we discussed in Sec. 4 in the scenario with two heavy neutral leptons.
A.1 Same sign leptons: We start with the decay rate for the process W + → + α + β q q mediated by two HNLs N 1,2 almost degenerate in mass, similar to that shown in Fig. 1, but without the initial quarks. The amplitude reads with p N = p W − p α . The channel with crossed α and β gives rise to the same amplitude, but it must be added incoherently to our process since the rate is dominated by on-shell N i , and thus the momentum of the first lepton is fixed by the 2-body decay kinematics. Therefore the two processes do not interfere and we can neglect for the moment the crossed channel. The only modification results in a factor of 2 in the rate. Defining U αN i = |U αN i | e iφ αi , δφ + = (φ α2 − φ α1 ) + (φ β2 − φ β1 ) and using the narrow width approximation (NWA), the squared matrix element becomes where ∆M 2 N ≡ M 2 N 2 − M 2 N 1 . Notice that, for the interference term, we used the NWA as follows: 1 which differs from the expression in Ref. [65], as discussed in Section 4. Assuming M N 1 M N 2 ≡ M N , Γ N 1 Γ N 2 ≡ Γ N and ∆M N ≡ M N 2 − M N 1 = 0, and considering that |U αN 1 | |U βN 1 | = |U αN 2 | |U βN 2 | ≡ |U αN | |U βN |, we get We observe that the squared amplitude in the case of only one sterile neutrino factorizes out.
Integrating over the phase space and after factorizing the decay width in the single HNL framework, it is straightforward to obtain Γ W + → + α + β q q = 2 1 + cos δφ + 1 1 + y 2 − sin δφ + y 1 + y where we have defined Notice that we obtain the same result of Ref. [63]. This is because we are using the NWA, which corresponds to consider on-shell HNL, leading to the same result after integrating over t, the time evolution of the intermediate HNL, from 0 to ∞. In fact the factor 1/Γ coming from the NWA is equivalent to the factor ∞ 0 e −Γt/2 2 = 1/Γ, and the interference part of Eq. (41) coincides with the finding of Ref. [63].
A.2 Different sign leptons: W + → + α − β q q In this appendix, the decay rate for the process W + → + α − β q q mediated by the two HNLs N 1,2 almost degenerate in mass is computed. The amplitude of this process is given by where, again, p N = p W − p α . With the help of FeynCalc [102], we obtain for the squared amplitude, with the K ij factors defined as Assuming as before that M N 1 ≈ M N 2 ≡ M N and Γ N 1 ≈ Γ N 2 ≡ Γ N , then K 11 = K 22 = K 12 ≡ K. Also, considering that |U αN 1 | |U βN 1 | = |U αN 2 | |U βN 2 | ≡ |U αN | |U βN |, we can simplify the expression notably with δφ − = (φ α2 − φ α1 ) − (φ β2 − φ β1 ). And using the narrow width approximation, we get Finally, it is straightforward to obtain the total decay width in terms of the decay width mediated by just one HNL, Γ W + → + α − β q q = 2 1 + cos δφ − 1 1 + y 2 − sin δφ − y 1 + y with y defined in Eq. (43).