Disambiguating Seesaw Models using Invariant Mass Variables at Hadron Colliders

We propose ways to distinguish between different mechanisms behind the collider signals of TeV-scale seesaw models for neutrino masses using kinematic endpoints of invariant mass variables. We particularly focus on two classes of such models widely discussed in literature: (i) Standard Model extended by the addition of singlet neutrinos and (ii) Left-Right Symmetric Models. Relevant scenarios involving the same"smoking-gun"collider signature of dilepton plus dijet with no missing transverse energy differ from one another by their event topology, resulting in distinctive relationships among the kinematic endpoints to be used for discerning them at hadron colliders. These kinematic endpoints are readily translated to the mass parameters of the on-shell particles through simple analytic expressions which can be used for measuring the masses of the new particles. A Monte Carlo simulation with detector effects is conducted to test the viability of the proposed strategy in a realistic environment. Finally, we discuss the future prospects of testing these scenarios at the $\sqrt s=14$ and 100 TeV hadron colliders.


Introduction
Despite the spectacular experimental progress in the past two decades in determining the neutrino oscillation parameters [1], the nature of new physics beyond the Standard Model (SM) responsible for nonzero neutrino masses is still unknown. A widely discussed paradigm for neutrino masses is the so-called type-I seesaw mechanism [2][3][4][5][6] which postulates the existence of heavy right-handed (RH) neutrinos with Majorana masses. The masses of the RH neutrinos, synonymous with the seesaw scale, is a priori unknown, and its determination would play a big role in vindicating the seesaw mechanism as the new physics responsible for neutrino mass generation. There are a variety of opinions as to where this scale could be [7], ranging from the left-handed (LH) neutrino mass scale of sub-eV all the way up to the grand unification theory (GUT) scale. The GUT approach, which is quite attractive, is mainly motivated by the straightforward embedding of the seesaw mechanism in SO(10) GUTs [3]. The seesaw scale in simple SO(10) models is near 10 14 GeV [7] making it quite hard to test in any foreseeable laboratory experiments. 1 There are also arguments based on naturalness of the Higgs mass which suggest the seesaw scale to be below 10 7 GeV or so [20,21]. It is therefore of interest to focus on the seesaw scale being in the TeV range since we have the Large Hadron Collider (LHC) searching for many kinds of TeV-scale new physics, which could also probe TeV-scale seesaw through the "smoking gun" lepton number violating (LNV) signal of same-sign dilepton plus dijet final states ± ± jj [22] and other related processes [23][24][25]. In addition, there are many kinds of complementary low energy searches for rare processes, such as neutrinoless double beta decay (0νββ) [26], lepton flavor violation (LFV) [27], anomalous Higgs decays [28][29][30][31] and so on which are sensitive to TeV-scale models of neutrino mass. For reviews on various phenomenological aspects of TeV-scale seesaw, see e.g., [32][33][34][35].
As far as the collider signals are concerned, in the case of SM seesaw, the Yukawa part of the model Lagrangian is given by where ψ L , H, N denote respectively the SM lepton and Higgs doublets and the RH neutrino singlet fields under SU (2) L , and C denotes the charge conjugation matrix. We have dropped the generation indices. 3 When H picks up a vacuum expectation value v to break the electroweak gauge symmetry, it leads to the usual type-I seesaw [2][3][4][5][6]. On the other hand, in the LRSM, based on the gauge group SU (2) L × SU (2) R × U (1) B−L , we have the additional gauge bosons W R , Z which play an important role in collider phenomenology, if W R mass is in the accessible range. The part of the Lagrangian relevant to our discussion is that involving the charged current and Yukawa interactions responsible for the seesaw matrix after SU (2) R × U (1) B−L symmetry is broken: There exist examples [40] where the small neutrino masses via type-I seesaw at TeV-scale can arise without excessive fine tuning of the LRSM parameters. 3 As in most of the collider studies of seesaw, we will assume later a single heavy neutrino giving the dominant contribution to the jj signal for simplicity.
where we have again omitted the flavor structure in the RH sector. Note that there is now a direct coupling of the RH neutrinos to the W R gauge boson which will distinguish its production mode at the LHC from that of the SM seesaw case. Assuming m N < m W R , as suggested from vacuum stability arguments [68], the basic collider signal arises from the on-shell production of the RH neutrinos accompanied by a charged lepton in the first stage, and the subsequent decay of N to jj final states (with = e, µ, τ , depending on the initial flavor of N ). The latter can go via the virtual intermediate state W R or physical on-shell W , assuming that m N > m W . 4 It is worth emphasizing that in the LRSM, the Majorana nature of the RH neutrinos inevitably leads to the LNV signature of ± ± jj [22], irrespective of the Dirac Yukawa couplings in Eq. (1.2). On the other hand, in the case of the SM seesaw, the strength of the dilepton signal crucially depends on the size of the heavy-light neutrino mixing parameters V N ∼ vh ν m −1 N , m N being the corresponding eigenvalues of the RH neutrino mass matrix M N in Eq. (1.1). In fact, the constraints of small neutrino masses from type-I seesaw have severe implications for the LNV ± ± jj signals [73][74][75]. Whether the dilepton signal can be of same-sign or mostly of opposite-sign depends on how degenerate the RH neutrinos are and to what extent they satisfy the coherence condition [76]. 5 Nevertheless, the general kinematic strategy presented in this paper is equally applicable to both same and opposite-sign dilepton signals. In this sense, its efficacy is not just limited to the type-I seesaw models, but also to many of its variants, such as the inverse [82,83], linear [8, 84,85] and generalized [17,86,87] seesaw models, which typically predict a dominant oppositesign dilepton signal. In fact, in light of the recent observations from both CMS [88] and ATLAS [89] indicating a paucity of ± ± jj events, the ± ∓ jj signal might turn out to be more relevant [17,18,90,91] in the potential discovery of parity restoration in future LHC data. Therefore, we will not specify the lepton charge in our subsequent discussion, and the signal will be simply referred to as the jj signal.
For m W R > m N > m W , there are four different sources in the LRSM for the origin of the jj signal at the LHC [47,92] (see Figure 1): where the first (LL) mode is the only one that arises in the SM seesaw via s-channel exchange of the SM W -boson (which is denoted by W L to justify the name LL), whereas all the four modes can arise in L-R models. These signals are uniquely suited to probe the Majorana and Dirac flavor structure of the neutrino seesaw and are therefore an important probe of the detailed nature of the seesaw mechanism. To this end, it is important to 4 For mN < mW , one can look for other interesting signals like displaced vertices at the LHC [69] or decay products of charged mesons [70][71][72], which we do not discuss here. 5 For instance, if the pattern of RH neutrino masses is hierarchical while maintaining large |V N | 2 [77] and satisfying the neutrino oscillation data, one can in principle have a ± ± jj signal [23][24][25][78][79][80][81]. The main point of this paper is to show that determining the kinematic endpoints of different invariant mass distributions, 6 e.g., m , m jj etc., provides a unique way to distinguish these scenarios, irrespective of the dynamical details. Note that given possible scenarios listed above, the invariant mass variables involving a single lepton suffer from the combinatorial ambiguity due to the unawareness of the order between the two leptons. A special prescription is taken for those variables. We then provide a systematic way of determining the kinematic endpoints of the modified variables resulting from such a prescription. Various ways of distinguishing the underlying new physics scenarios by making use of the sharp features of kinematic variables have been studied in several recent works, e.g., in [94][95][96][97] in the context of dark matter stabilization symmetries. In particular, kinematic endpoints are robust against detailed dynamics of the underlying physics such as spin correlation, i.e., specific and extreme kinematic configurations are relevant to the kinematic endpoints without the need to know the full dynamical details of the process. Moreover, they are typically protected from event selection cuts unless such cuts are customized to reject the events corresponding to the kinematic endpoints. These observations motivate us to utilize the endpoints of various invariant mass variables to distinguish between the seesaw signals (1.3)-(1.6). We emphasize that different scenarios typically give rise to different event topologies, resulting in the kinematic endpoints having distinctive dependencies upon underlying mass parameters. The identification of specific relationships among them will enable us to uncover the relevant model parameter space and may eventually lead to a measurement of the masses of the involved heavy particles. Hence, we argue that this method is more robust than the dynamical variables like spin or angular correlations, as proposed earlier to distinguish the seesaw signals [47,98].
The rest of the paper is organized as follows. We first define our notations to be used in subsequent sections, followed by a discussion about the general strategy, in Section 2. We then give detailed derivations for the kinematic endpoints of various invariant mass variables case-by-case in Section 3. In Section 4, we discuss various ways for the topology disambiguation and mass measurements of the new particles. In Section 5, we illustrate the utility of the endpoint technique for the L-R seesaw signals by performing a Monte Carlo simulation including detector effects. In Section 6, we discuss the L-R seesaw phase diagram for collider studies and the future prospects of probing the LRSM parameter space at the √ s = 14 TeV LHC, as well as its complementarity with other low-energy probes. We also derive the first LRSM sensitivity contours for the planned √ s = 100 TeV pp collider. Our conclusions are given in Section 7. In Appendix A, we illustrate some invariant mass distributions for the LRSM.

Notations and general strategy
To accommodate all possibilities for an jj final state in a model-independent manner, we generalize the scenarios introduced in Eqs. (1.3)-(1.6) according to their decay topologies in conjunction with symbolic notations. We begin with a three-step cascade decay sequence of a heavy resonance C: where n and f are correspondingly identified as "near" and "far"-lepton with respect to the particle C to identify their relative location to each other. The associated event topology is explicitly diagrammed in Figure 2(i), being henceforth denoted as Case (i). If particle A is heavier than particle B, then the latter directly decays into f jj via a 3-body decay as shown in Figure 2(ii), which is denoted as Case (ii). In an analogous manner, one can imagine the situation where particle B is heavier than particle C so that the latter decays into two leptons and particle A via a 3-body decay. The relevant diagram is shown in Figure 2(iii) and denoted as Case (iii). Another possibility is the situation where particle C directly decays into two leptons and two jets via a 4-body decay, as shown in Figure 2(iv) Topology Region(s) Model Scenario(s) Table 1. Relations of the event topologies (i)-(vii) shown in Figure 2 with the relevant kinematic regions labelled as R 1 − R 7 , as well as their mapping with the possible scenarios in the LRSM. The subregions in parenthesis, e.g., R 1(1) and R 1 (2) , mean that they correspond to the same event topology but with different kinematic endpoints, as explained later in Section 3 [cf. Eqs. and labelled as Case (iv). Finally, we consider the situations where particle C is off-shell (denoted as C * ) for Cases (i) through (iv), and their counterparts are respectively exhibited in Figure 2(v) to 2(viii), and labelled as Cases (v) to (viii). Note that Case (viii) is a very unlikely scenario, as C is always assumed to be heavier than the jj system, so we do not consider it any further. Provided with the L-R seesaw models described in the previous section, one can easily correlate these event topologies with various scenarios arising therein by identifying N as particle B and various combinations of gauge bosons (W L , W L ), (W R , W R ), (W R , W L ), and (W L , W R ) as particles (C, A), respectively, depending on the underlying scenario. To accommodate all possible mass hierarchies, we introduce a subscript h (l) meaning that N is heavier (lighter) than the subscripted gauge boson, e.g., R l L h means m W R > m N > m W . For LL, RR, RL, and LR, each letter accepts either l or h subscript so that (naively) 16 scenarios would be possible. However, the subscripts in four scenarios such as L l L h , L h L l , R l R h , and R h R l imply a contradictory mass spectrum, i.e., the same gauge boson (either W L or W R ) cannot be heavier and lighter than N simultaneously. In addition, we have only assumed m W R > m W , as required to satisfy the experimental constraints from direct searches at the LHC [88,89], as well as the indirect constraints from K L − K S mass difference [99][100][101]. Therefore, the scenarios implying m W > m W R , i.e., R h L l and L l R h , are not allowed. As a result, only 10 different combinations are possible. In principle, the particle identified as C can be forced to be on-shell (assuming that it is within the kinematic reach of the collider experiment), so that all 10 possibilities can be assigned to the first four event topologies, listed in Table 1. Depending on the underlying model details, particle C can also be off-shell, if m C < m B , or even N can be off-shell. The relevant identifications to the last three topologies are also tabulated in Table 1. Note that the scenarios listed in Eqs. (1.3)-(1.6) can be obtained readily by imposing the additional assumption that m W R > m N > m W , as shown by the colored text in Table 1. A detailed topology identification of these scenarios with respect to the kinematic regions R 1 to R 7 will be discussed in Sections 4 and 5.
With the final state of n f jj, one can come up with eight non-trivial invariant mass variables, namely, m , m nj , m f j , m jj , m j , m njj , m f jj , and m jj , where we omit the subscript on the lepton for the variables involving both leptons. We remark that the two jets are emitted from the same particle for all eight cases, and thus we do not distinguish one from the other. An immediate experimental challenge is that one does not know a priori which lepton comes first with respect to the C particle, although we theoretically label them as n and f . This innate combinatorial issue particularly becomes a major hurdle to form invariant mass variables involving a single lepton. For example, m nj or m f j are not valid experimental observables again because n and f are not discernible event-by-event.
Given an invariant mass variable having such a combinatorial ambiguity, a possible prescription to resolve it is to evaluate both invariant mass values for a fixed set of jets and divide them into the bigger and the smaller ones, labelled by the superscripts ">" and "<", respectively. Taking the example of the invariant masses formed by a lepton and a jet, we compute m nj and m f j for a given j and assign the bigger (smaller) one to m > j (m < j ). A similar rule is applied to the invariant mass variables constructed by a lepton and two jets, yielding m > jj and m < jj . With this prescription, we end up with the following eight experimentally valid observables: We basically measure the kinematic endpoints of these eight variables to identify the underlying decay topology by examining their interrelationship and kinematic features which are the main subjects of the next section. 7

Derivation of kinematic endpoints
In this section we derive the analytic expressions for the kinematic endpoints of the invariant mass variables listed in Eqs. (2.2)-(2.4). We first elaborate the detailed steps in Figure 3. The event topology for Case (i) (left panel) and the relevant kinematics described in the rest frame of particle B (right panel). Here, α is the polar angle of f with respect to the direction of n (or equivalently, particle C) while θ is the polar angle of j (black arrow) with respect to the direction of B in the rest frame of particle A. β and φ are azimuthal angles for n and j (black arrows) around the horizontal axis defined by f .
obtaining the final expressions when particle C is on-shell, and later we briefly mention how the relevant results can be applicable to the case of off-shell C. Our main focus is the kinematic endpoints that are insensitive to the details of relevant decays, and thus we do not hypothesize any specific matrix element, i.e., we deal with only the phase-space structure for the associated decay kinematics.
This scenario can be represented by a three-step on-shell cascade decay of a heavy resonance C, as shown in Figure 3(a). We emphasize that it is rather convenient to describe the associated kinematic configuration in the rest frame of particle B, as shown in Figure 3(b). 8 Naively, the total degrees of freedom are four, but one angle β can be dropped by taking into account the azimuthal symmetry of the system. The mass hierarchy in this case determines the baseline inequalities defining the region of interest: where R ij denotes the mass ratio between the massive states i and j: R ij ≡ m 2 i /m 2 j .

2-body invariant mass variables
We begin with the invariant mass of two leptons m , for which the kinematic endpoint is well-known: Although angle θ is defined with respect to one jet of the two, one could develop the parallel argument with respect to the other jet. Since θ is a quantity measured in the rest frame of particle A, exactly the same argument goes through with cos θ → − cos θ. However, this does not provide additional information due to the fact that cos θ (− cos θ) spans +1 (−1) to −1 (+1) as θ increases, and in turn, any similar prescription which would be applied to jet-induced combinatorics for a fixed set of leptons does not enable us to obtain any further independent information.
Trivially, m 2 jj is given as a resonance of m 2 A : In order to derive the expressions for m >,max j and m <,max j , we first rewrite the relevant two invariant masses in terms of angular variables introduced in Figure 3(b): with η, x m , and z m defined as follows: Note that the invariant mass variables in Eqs. (3.4) and (3.7) are normalized by m 2 C for later convenience.
Since the angular variables α and φ are irrelevant to variable x, we first maximize the variable z over them to end up with an expression of z as a function of only cos θ: which is nothing but a straight line in cos θ with a negative slope. To obtain the maxima for m > j and m < j , we compare Eqs. (3.4) and (3.7) while varying cos θ. Since x| cos θ=1 = 0 while z| cos θ=1 = 0, two topologically different cases arise as shown in Figure 4. Each line describes the maximum x or z value for a given cos θ. For x m ≥ z m (left panel), one can identify which invariant mass is greater or smaller for any given event, i.e., for a fixed cos θ. Such a comparison can be conducted for all values of cos θ, and the trajectories for the greater and the smaller are explicitly delineated by blue-solid and red-dashed arrows respectively. The kinematic endpoint can be obtained simply by reading off the maximum value of each line. For the current example (left panel), they are x m and z m , and vice versa for x m < z m (right panel). Thus, the invariant masses can be formulated as follows: with R AB ∈ (0, 1]. We denote the first and the second regions in Eq. (3.8) as R 1(1) and R 1 (2) , respectively, and their link to specific model scenarios are tabulated in Table 1.

3-body invariant mass variables
For the kinematic endpoint for m 2 j , we simply use the relevant expression from Ref. [103]: Considering the other two variables m > jj and m < jj , we remind the reader that the two jets originate from a common mother particle A, and thus they are equivalent to m > A and m < A , i.e., they are reduced to 2-body invariant mass variables involving a massive visible particle. To obtain their analytic expressions, we first find the formulae for m 2 f A and m 2 nA . The former is nothing but the resonance B, i.e., m 2 f A = m 2 C R BC . The maximum of the latter can be easily obtained by going through the steps leading to (m max ) 2 [cf. Eq. (3.2)]: Since m 2 f A is fixed whatever m 2 nA is given, the kinematic endpoints for the variables of interest can be obtained, as follows: where the first and the second regions correspond to R 1(1) and R 1(2) , respectively.

4-body invariant mass variable
The 4-body invariant mass variable m jj is trivially given by the mass of resonance C: Figure 5. The event topology of Case (ii) (left panel) and the relevant kinematics described in the rest frame of particle B. Here α is the polar angle of n with respect to the direction of the composite system of f and j (black arrow) while θ is the polar angle of particle f with respect to the direction of the composite system in the rest frame of particle B. β and φ are azimuthal angles for particles n and f about the axis extended by j (red arrow) and the composite system.
Applying the formulae to Case (v) Once the particle C is off-shell, as in Case (v) shown in Figure 2, one can interpret that its mass is given by the center of mass energy √ŝ . Therefore, it is possible to reuse all analytic expressions thus far by replacing m 2 C byŝ max = s. The chance of reachingŝ max is, however, so small that any associated endpoints depending on s are not typically saturated. Therefore, only the variables having no dependence on s are reliable: m 2 jj , one of m 2 j 's, and one of m 2 jj 's. Since s is assumed much larger than the mass of particle B, R BC → 0, 9 so that the corresponding region (denoted by R 5 in Table 1) is defined by R BC = 0 and 0 < R AB < 1 . (3.12) This scenario can be represented by a two-step cascade decay of a heavy resonance C, where the second step is proceeded via a 3-body decay. The corresponding diagram and kinematic configuration are delineated in Figure 5. It turns out that it is convenient to do the analysis in the rest frame of particle B as shown in Figure 5(b). Just like Case (i), the angle β is irrelevant to the description of the system due to the azimuthal symmetry. The mass hierarchy for this case defines the baseline inequalities determining the corresponding region as

2-body invariant mass variables
Since particle B decays into three final states via off-shell A, m 2 jj is given by a distribution, not a resonance as in Case (i), and the analytic expression for its endpoint is 14) The kinematic endpoint for m 2 can be found in Ref. [104], i.e., When it comes to the analytic expressions for m >,max j and m <,max j , we go through a similar argument to that in Case (i).
We observe that x m is simply given by R BC . Unlike the previous case, cos θ is irrelevant to x so that one can maximize z even over cos θ as well as φ and α, yielding z = z m which is nothing but a horizontal line in x . Clearly, only two cases are available, and the respective expressions are and we denote the first and the second regions as R 2(1) and R 2(2) , respectively.

3-body invariant mass variables
To obtain the expression for m max j 2 , the result in Ref. [104] can be reused, simply yielding For the other two variables, we again follow a similar argument to the previous case: Applying the formulae to Case (vi) As in Case (i), we can simply reuse all expressions derived thus far along with replacement of m 2 C by s for an off-shell C. Again, the chance of havingŝ max = s is so small that the kinematic endpoints depending on s are not typically saturated. Therefore, only the variables having no dependence on s are reliable: m 2 jj , one of m 2 j 's, and m < jj 2 . As in the previous case, s is assumed to be much larger than the mass of particle B, and thus the inequalities defining the corresponding region (denoted by R 6 ) are given by Figure 6. The event topology of Case (iii) (left panel) and the relevant kinematics described in the rest frame of particle C. α is the polar angle of n with respect to the direction of the composite system of n and f while θ is the polar angle of j (black arrow) with respect to the direction of such a composite system at the rest frame of particle C. β and φ are azimuthal angles for n and j (black arrow) about the axis extended by j (red arrow) and the composite system.

Case
This scenario can be represented by a two-step cascade decay of a heavy resonance C, where the first step is proceeded via a 3-body decay. The corresponding diagram and kinematic configuration are delineated in Figure 6. It turns out that it is convenient to do analysis in the rest frame of particle C, as shown in Figure 6(b). 10 Just like the previous two cases, azimuthal angle β is irrelevant to the description of the system due to the azimuthal symmetry. Though particle B is heavy enough to be integrated out, R AC should be still less than 1, i.e., R AC = R AB R BC < 1. Therefore, the baseline inequalities defining the corresponding region (denoted by R 3 ) are (3.20)

2-body invariant mass variables
Since particle A undergoes a 2-body decay, m 2 jj is simply given by the resonance of particle A, i.e., m 2 jj = m 2 C R AC . Due to the 3-body decay of C, the endpoint of m 2 is given by For the purpose of deriving the analytic expressions of kinematic endpoints for m > j and m < j , we first express the momenta for all visible particles in the rest frame of particle C: p µ j = m A 2 (cosh η 2 + sinh η 2 cos θ, sinh η 2 + cosh η 2 cos θ, sin θ cos φ, sin θ sin φ) , (3.24) where η 1 and η 2 are defined as With these definitions, one can prove the following useful relations: where λ denotes the kinematic triangular function defined as λ(x, y, z) ≡ x 2 + y 2 + z 2 − 2(xy + yz + zx) . (3.27) The symmetry of n ↔ f enables us to have 28) which is maximized at cos θ = −1 and y = 0, i.e., From this constraint, we find that m > j can be maximized when either x or z vanishes, and m < j can be maximized when x = z. Therefore, we have (3.30)

3-body invariant mass variables
To obtain the kinematic endpoint for m 2 j the relevant result in Ref. [105] can be reused: Note that the second expression corresponds to the kinematic configuration where the two visible particles are moving in the same direction with equal energy while particle A is moving in the opposite direction.
Applying the formulae to Case (vii) As in the previous two cases, we again reuse all expressions derived thus far with the simple replacement of m 2 C by s. In this case, the only variable having no dependence on s is m jj and all other variables unreliable as the chance of havingŝ max = s is so small that the kinematic endpoints depending on s are not typically saturated. With the assumption that s is much larger than any of on-shell particle masses, we have R AC → 0. As R BC ≥ 1, R AB should be close to 0, giving the inequalities defining the corresponding region (denoted by R 7 ) as follows: R BC ≥ 1 and R AB = 0 . (3.33)

Case (iv): m A , m B ≥ m C
This scenario can be represented by a single 4-body decay of a heavy resonance C. A convenient frame for the relevant analysis is at the rest frame of particle C. The mass hierarchy of this case enables us to have the following inequality defining the corresponding region (denoted by R 4 ):

2-body invariant mass variables
Due to the 4-body decay feature of this process, m 2 jj and m 2 have an identical distribution whose kinematic endpoint is given by For the maximum of m > j , we imagine the situation where one of two 's and one of two j's are so soft in the rest frame of particle C that all mass-energy of C is split into the other and the other j. Similarly, for the maximum of m < j , we imagine the situation where one of j is extremely soft, the two leptons are emitted in the same direction, and the other j is emitted in the opposite direction so that m 2 nj = m 2 f j and m 2 = 0. Therefore, we have

3-body invariant mass variables
The expression for the m 2 j kinematic endpoint is simply given by

Summary
Let us collect the formulae we have derived so far for the kinematic endpoints and rearrange them with respect to invariant mass variables in Eqs. (2.2)-(2.4). To make a connection of the formulae with specific LRSM scenarios, we again refer the readers to Table 1.
where the expressions in the square brackets indicate that the associated distributions appear as a resonance peak. Since the mass parameter m 2 C merely determines the overall scale, the relevant parameter space can be divided in the  Table 1. Additional subscripts in the parenthesis denote subdivisions in the associated topology. as illustrated in Figure 7. As introduced earlier, R 1 -R 7 correspond to different decay topologies diagrammed in Figure 2 and tabulated in Table 1 [see also Eqs.  .17)]. Note that regions R 5 , R 6 , and R 7 are represented by one-dimensional strips unlike the others. This is because the basic assumption that the center of mass energy s is much larger than any of the masses of on-shell particles allows either R AB or R BC to vanish [cf. Eqs (3.12), (3.19), and (3.33)].

Topology disambiguation and mass measurement
In the busy environment of a real-life hadron collider experiment, leptons are much betterreconstructed than jets. Thus, from the precision measurement aspects, the best among the eight invariant mass variables listed in Eqs. (2.2)-(2.4) is m . It therefore follows that as a minimal choice, the dilepton invariant mass variable can be used for topology disambiguation. For example, Ref. [47] examined the shape of the dileption invariant mass distribution for the decay topologies (1.3)-(1.5). This was based on the observation that different decay topologies involve different spin correlations [98] which affect the m spectrum. However, typical challenges in shape analysis are the facts that (relatively) large statistics is required and more importantly, the relevant shape would be severely distorted by the detailed dynamics and the predicament of hard cuts essential to suppress the relevant SM background. 11 On the contrary, the kinematic endpoint does not suffer from not only the details of the associated decay but the cuts to suppress backgrounds because it depends only on some specific kinematic configurations. Also, the need for "local" information near the endpoint typically demands less statistics than the shape analysis. In any case, we lose useful information encoded in the shape of the distribution, and the measurement of m endpoint by itself does not enable us to distinguish a certain topology from others. Hence, we need to supplement this with other observables.

Constraints and Correlation matrix
Following the basic idea of involving as few jets as possible in the observables to be used, the next best ones are the invariant mass variables having only a single jet, i.e., m > j , m < j , and m j . However, including m , this makes the number of observables greater than the number of unknown mass parameters (m C , m B , and m A ), namely, the system is now over-constrained. 12 We find several characteristic sum rules/constraints among the four invariant mass variables mentioned above, which are listed below in terms of the notations introduced in Eqs. C6: Only c is a well-defined endpoint for R 5 , R 6 , and R 7 .
C7: For R 1(2) and R 5 , m > j has a "foot" structure in it, i.e., a cusp is developed in the middle of the distribution. The cusp position should be the same as c.
We also provide a correlation matrix between regions and sum rules/constraints given above in Table 2. Here the symbol " √ " implies that the relevant sum rule/constraint should be obeyed for the region of interest. We see that in principle, it is possible to discern all regions but R 6 and R 7 using C1-C7.
Of course, more sum rules/constraints can be utilized, once we allow the invariant mass variables involving more than 1 jet, which enable us to distinguish even R 6 and R 7 . Instead of exhausting all possibilities, we enumerate some of them below: C8: d appears as a resonance peak for R 1(1) , R 1(2) , R 3 , R 5 , and R 7 .
11 Despite such potential difficulties, the invariant mass shape can in some cases be sufficient to completely determine the new particle mass spectrum, including the overall mass scale [96]. 12 For some specific cases of LRSM, either A or C, or both can be the SM W boson. One would then expect the number of unknown mass parameters reduced to two or one. However, since we do not know which underlying scenario governs the channel of our current interest, we generically treat the masses of A, B, C as unknown parameters in our discussion. Table 2. Correlation matrix between regions and sum rules/constraints. " √ " implies that the relevant sum rule/constraint should be satisfied for the region of interest. C1-C7 involve only the invariant mass variables having 0 or 1 jet, while the others involve ≥ 2 jets.
It may be noted here that Refs. [106,107] studied a way of distinguishing event topologies (i) and (ii) in Figure 2 by the existence of resonance peaks, which is relevant to C8 and C9.

Mass measurement of new particles
Once different regions are determined with the methods explicated in the previous section, the mass measurement of on-shell particles, or equivalently, the measurement of R ij 's together with mass measurement of the first on-shell particle mass can be readily performed in terms of kinematic endpoints of various invariant mass distributions. We first provide the inverse formulae using only a, b, and c for various regions but R 5 and R 7 .
(4.5) For R 5 and R 7 , the relevant mass parameters can be extracted by using other observables involving more than one jet. For example, We should remark that Eqs. (4.1)-(4.9) are not unique ways of writing inverse formulae. For most of the regions, the interrelationships among the eight available kinematic endpoints are over-constrained due to less number of unknown mass parameters than the employed variables. Although the above-listed inverse formulae are expressed with a certain set of kinematic endpoints, others can as well be utilized for cross-checks.

Application to the L-R seesaw
In this section, we apply the general kinematic endpoint technique discussed in the preceding section to the specific case of L-R seesaw. Our aim is to illustrate the distinction between the invariant mass distributions for the scenarios mentioned in Eqs. In terms of the model scenarios listed in Table 1, they correspond to L h L h , R l R l , and R l L h , respectively. One can therefore correspondingly relate them to event topologies (v), (ii), and (i), and hence, to regions R 5 , R 2 and R 1 , respectively in Figure 7. Parton-level events are generated with MadGraph aMC@NLO [108] and the parton distribution functions (PDFs) of protons are evaluated by the default NNPDF2.3 [109]. To describe parton showering and hadronization, the events are streamlined to Pythia6.4 [110]. The relevant output is subsequently fed into Delphes3 [111] interfaced with FastJet [112] for describing detector effects and finding jets. All the simulations are conducted for a √ s = 14 TeV pp collider at the leading order. In the first two benchmark scenarios, jets are formed using the anti-k t algorithm [113] with a radius parameter R = 0.4. For the last benchmark scenario, we point out that the on-shell W bosons tend to be highly boosted due to a large mass gap between N and W . The large boost of W leads to a high collimation of the jets from its decay, thus requiring jet substructure techniques to resolve the subjets in each two-prong "W "-jet. To this end, jets are initially clustered by the Cambridge-Aachen algorithm [114,115] with a jet radius R = 1.2, and the resulting "fat" jets are further processed by the Mass Drop Tagger method [116].
We emphasize that the main purpose of this simulation study is to see whether or not the proposed endpoint strategies are viable in the presence of detector effects. In this sense, including any potential backgrounds is beyond the scope of this paper, so we simply consider the signal component. Also, the precise measurement of kinematic endpoints typically requires large statistics; therefore, an arbitrarily large number of signal events are generated to minimize any statistical fluctuation. Finally, event selection is executed with a minimal set of selection criteria for simplicity. Of course, in the presence of backgrounds, one would impose a more stringent set of cuts to suppress them, but considering the fact that kinematic endpoints are typically least affected by cuts, we expect that the endpoint behaviors in our scheme will be similar to those under more sophisticated selection criteria.
Particle acceptance/isolation and detector geometry are basically performed according to the default parametrization in Delphes3 [111]. We then choose the events having exactly two leptons (of either e or µ flavor) and ≥ 2 jets in the final state. For the second scenario (BS2), the first two hardest jets are used to construct relevant invariant masses, while for the other two scenarios (BS1 and BS3), the W mass window is employed because the relevant jets were emitted from an on-shell W boson. More specifically, we evaluate the dijet invariant mass for all possible combinations among observed jets and choose two jets satisfying a slightly tight condition of |m jj − m W | < 12 GeV. If there exist multiple combinations, then we simply take the one with the smallest difference. Although some events can contain more than 2 reconstructed jets, we evaluate the jet-involving invariant mass values using the two jets identified in the above-explained way, so that in every single event we have the same number of invariant mass values.
We show various 2-and 3-body invariant mass distributions for the three benchmark scenarios in Figure 8. 13 Top, middle, and bottom panels correspond to BS1, BS2, and BS3, respectively. In each row of panels, the left panel shows the 2-body invariant mass distributions m (red dashed), m > j (blue dot-dashed), and m < j (green solid) while the right panel shows the 3-body invariant mass distributions m j (red dashed), m > jj (blue dot-dashed), and m < jj (green solid). Black dashed lines indicate the relevant theory predictions for the kinematic endpoints derived in Section 3. Note that we provide unit-normalized distributions for illustration since our focus is the structure of kinematic endpoints. Based on the analysis scheme described above, one can easily infer the corresponding relative weights. We observe that for BS1 and BS2 the invariant mass variables involving ≤ 1 jet, i.e., m , m > j , m < j , and m j , are reasonably well-matched to the associated theoretical predictions in spite of detector effects such as jet energy resolution, smearing and contamination from initial and final state radiations (ISR/FSR). 14 In contrast, once more jets are involved, 13 Dijet invariant mass distributions are not displayed in the figure, as it develops a trivial distribution, i.e., a resonance peak, for BS1 and BS3.
14 Note that for BS1 the kinematic endpoints for m , m > j , and m j are ill-defined because their values can be arbitrarily large up to the associated kinematic limit determined by the available center of mass energy. i.e., m > jj and m < jj , kinematic endpoints are more smeared so that identifying their correct positions would be rather challenging. When it comes to BS3, the situation becomes more challenging, and even for 2-body invariant mass distributions, the endpoints either suffer from relatively large smearing (green and red histograms) or are unsaturated (blue histogram). In particular, the information of subjets in the W -jet is typically less accurate than that of regular jets, thus rendering the distributions more smeared. Certainly, these phenomena can be improved by better jet energy measurement and ISR/FSR identification (see, for example, Ref. [117] studying the impact of various detector effects upon the distributions of mass variables). Moreover, the boosted jet technique is a research field that is being actively investigated and developed. We therefore expect that locating the true kinematic endpoints will be advanced in the future, which will help us in determining the new particle masses more accurately.

L-R seesaw phase diagram
In this section, we present some future prospects of probing the L-R seesaw parameter space using the different collider signals mentioned earlier. To be concrete, we only focus on the mass hierarchy m W R > m N > m W which can be effectively probed at the LHC via the jj final state. In this case, as mentioned in Eqs. (1.3)-(1.6), there are four classes of Feynman diagrams for the jj signal ( Figure 1) [47]. The LL channel is a clear probe of the seesaw matrix in both SM seesaw and L-R seesaw models. However, its effectiveness solely relies on the heavy-light neutrino mixing parameter |V N | 2 , and is limited to heavy neutrino mass M N only up to a few hundred GeV [23][24][25][78][79][80][81]. Experimentally, the mass range M N = 100-500 GeV has been explored at the √ s = 8 TeV LHC for = e, µ [89,118], and direct upper limits on |V N | 2 of the order of 10 −2 -10 −1 have been set. We note here that the current indirect limits from a global fit of electroweak precision data (EWPD), lepton flavor violation, and lepton universality constraints [119] are roughly one to two orders of magnitude stronger. For a bird's-eye view of other complementary limits and the future sensitivities, see e.g., [34].
In the RR channel [22], the RH neutrino is produced on-shell due to its gauge interaction with W R [cf. Eq. (1.2)] and subsequently decays into a 3-body final state via an off-shell W R . This diagram only relies on the gauge coupling g R and is independent of V N ; therefore, it gives the dominant contribution for small V N which is of course the naive expectation in the "vanilla" type-I seesaw case. Using this channel, LHC exclusion limits are derived in the (m N , m W R ) plane, and currently exclude m W R up to 3 TeV assuming the equality between the SU (2) R and SU (2) L gauge couplings, i.e., g R = g L [88,89]. 15 In general, the signal cross section will be scaled by a factor of (g R /g L ) 4 , and hence, the corresponding limit on m W R could be weaker if g R < g L . In any case, being independent of the Dirac neutrino Yukawa coupling, the RR process does not probe the complete seesaw matrix.
The RL and LR contributions, on the other hand, necessarily involve the heavy-light neutrino mixing. In fact, the RL diagram could give the dominant contribution to the jj signal, if the mixing |V N | 2 is not negligible and/or the W R gauge boson is not too heavy [47]. There are two reasons for this dominance: (i) it leads to a production rate σ(pp → W R → N ± ) which is independent of mixing and only suppressed by (g R /g L ) 4 (m W /m W R ) 4 (as in the RR case), and can therefore dominate over the LL contribution which depends on |V N | 2 g 4 L ; (ii) the RH neutrino in this case has a 2-body decay: N → ± W L → ± jj (as in the LL case) which is not suppressed by the phase space, unlike in the RR case with a 3-body decay of N . Hence, for a sizable range of the mixing and RH gauge boson mass, the RL mode is expected to be dominant for the jj signal at the LHC and could constitute 15 Similar lower limits on mW R are also obtained from low-energy flavor changing neutral current observables [102]. a clear probe of the full seesaw matrix. There exist classes of L-R seesaw models where such large mixing parameters can be realized without much fine-tuning [40]. The remaining possibility, namely, the LR contribution is doubly suppressed by the heavy-light mixing as well as phase space, and hence, always smaller than one of the other three contributions discussed above. So this is not relevant for the experimental searches of L-R seesaw.
The regions of dominance for the different contributions discussed above are shown in Figure 9 by various shaded regions (blue for LL, green for RL and red for RR) as a function of the mixing parameter in the electron sector for a fixed value of the heavy neutrino mass m N = 1 TeV and assuming g R = g L . The vertical (brown) solid line in Figure 9 shows the 95% CL direct limit on m W R from the √ s = 8 TeV LHC [118], whereas the vertical dashed line shows a projected lower limit from √ s = 14 TeV LHC with 300 fb −1 integrated luminosity [41]. For comparison, we also show the current 90% CL upper limit on the mixing parameter from a recent global fit to the EWPD [119] (horizontal dotted line). All the signal cross sections for LL, RL and RR used in Figure 9  For simplicity, we assume all the non-standard Higgs bosons in the LRSM to be heavier than W R , so that the total width of W R is mostly governed by the masses of W R and N .  Figure 10. L-R seesaw phase diagram for the jj signal at the LHC with |V eN | 2 = 10 −6 . The labels are the same as in Figure 9. In addition, we have shown the exclusion region due to LFV constraints from MEG (solid orange) and the future limit from MEG2 (dashed orange). The grey shaded region (upper left corner) corresponds to m N > m W R which cannot be probed by the jj signal at the LHC, but accessible to the low-energy searches.
The presence of Majorana neutrinos in the LRSM also gives several additional contributions to the low-energy LNV process of 0νββ [3,40,[50][51][52][53][54][55][56][57][58][59][60][61][62]. In the limit of large m W R , the only dominant contributions are due to LH current exchange, which gives an upper limit on the mixing parameter, independent of m W R , as shown by the solid magenta line in Figure 9. On the other hand, for very light RH gauge bosons, the purely RH current exchange diagram becomes dominant over the rest and puts a lower bound on m W R for a given m N from the non-observation of 0νββ, irrespective of the value of the mixing parameter, as shown by the vertical portion of the solid magenta curve. Here we have used the 76 Ge isotope and the corresponding combined 90% CL limit on the 0νββ half-life from GERDA phase-I: T 0ν 1/2 > 3 × 10 25 yr [120] for illustration. For the nuclear matrix elements, we use the maximum values from a recent SRQRPA calculation [121], so as to obtain the minimum half-life predictions. The dashed magenta curve shows the future sensitivity of multi-ton scale 0νββ experiments (LSGe) with 76 Ge isotope, such as the proposed Majorana+GERDA experiment with an ultimate limit of T 0ν 1/2 > 10 28 yr [122]. The phase diagram shown in Figure 9 can be easily translated to the more familiar (m W R , m N ) parameter space, as shown in Figure 10. The current 95% CL exclusion contour from √ s = 8 TeV LHC [118] (LHC8) and the future sensitivity of √ s = 14 TeV LHC with 300 fb −1 luminosity (LHC14) [41] are also shown. Here we have fixed the mixing parameter |V eN | 2 = 10 −6 for illustration. Due to the relatively small value of mixing, the LL channel is no longer relevant, and either RL or RR channel is the dominant one in the entire L-R seesaw parameter space of interest. In particular, for smaller m N values, the RR channel becomes less efficient, compared to the RL channel, due to the phase space suppression factor of m 5 N /m 4 W R in the 3-body decay rate Γ(N → jj). Increasing the value of the mixing parameter will further enlarge the RL dominance region. Thus, a combination of RR and RL modes provides a better probe of the L-R seesaw, as compared to the RR mode alone. It is worth noting here that the 0νββ searches provide a complementary probe of the L-R seesaw parameter space, especially for the m N > m W R regime, which is kinematically inaccessible in the jj channel at the LHC. We should note here that the 0νββ analogs of the LL, RL and RR modes will give rise to different angular distributions, which can in principle be measured in the proposed SuperNEMO experiment [123]. Similarly, one can use polarized beams in a linear collider to distinguish between the different contributions in L-R seesaw [124]. These are complementary to the kinematic endpoint method proposed here for a hadron collider.
Another complementary low-energy probe of the LRSM is through the LFV processes [45,56,62,66,67]. In particular, the µ → eγ decay rate gets an additional contribution from the purely RH current and is currently constrained to be Br(µ → eγ) < 5.7×10 −13 at 90% CL by the MEG experiment [125]. Assuming a maximal mixing between the RH electron and muon sectors of the heavy neutrino mass matrix, we translate this limit into an exclusion region in the (m W R , m N ) parameter space, as shown by the shaded region above the solid orange curve in Figure 10. The projected limit of Br(µ → eγ) < 10 −14 from the future upgrade of MEG experiment [126] could probe most of the remaining parameter space shown in Figure 10. This clearly illustrates the importance of a synergistic approach at both energy and intensity frontiers in testing the L-R seesaw paradigm in future.
Since the high energy physics community has seriously started thinking about a nextgeneration √ s = 100 TeV hadron collider [127], we feel motivated to present the sensitivity reach of such a machine in the context of L-R seesaw. As before, we generate the signal cross sections at the parton level using MadGraph aMC@NLO [108] with the default NNPDF2.3 [109] PDF sets. We have used the following conservative selection cuts in our event simulation: Our results for the projected 2σ exclusion region with 1 ab −1 integrated luminosity at √ s = 100 TeV are shown in Figure 11 (blue dotted curve), which suggests that one could probe RH gauge boson masses up to 32 TeV, in good agreement with previous studies [48]. This sensitivity reach extends to regions of the LRSM parameter space not accessible by even future low-energy searches at the intensity frontier, as demonstrated by a comparison with the MEG-2 projected limit (orange dashed curve) in Figure 11. Moreover, one can access even smaller mixing parameters at the collider through the RL mode, as illustrated here for |V N | 2 = 10 −8 which still gives a dominant RL contribution for relatively smaller m N values. Note that the shape of our exclusion contour for the √ s = 100 TeV case is not exactly similar to the √ s = 8 or 14 TeV contours in the low m N region, mainly because we have not applied any specialized selection cuts on the invariant masses and have taken the tagging efficiencies to be 100%, which is usually not the case in a realistic hadron collider  Figure 11. Sensitivity reach of a futuristic √ s = 100 TeV collider for probing the L-R seesaw parameter space through the jj signal. Other labels are the same as in Figure 10. environment. Nevertheless, our parton-level estimates should serve as a rough guideline for more sophisticated studies in future.

Conclusion
In summary, we have pointed out a clean and robust way to distinguish between different mechanisms for the production of dilepton plus dijet final states at a hadron collider using the kinematic endpoints of various invariant mass distributions. We derived analytic expressions for these kinematic endpoints, under a minimal set of assumptions : (i) no invisible particles are involved in the relevant process (i.e., no missing transverse momentum at the parton level), (ii) all the final state particles originate from a common mother particle, and (iii) the two jets are the decay products of the same particle. We then provided various criteria to distinguish the possible scenarios yielding the same jj collider signature from one another. As a "spin-off", the potential determination of the masses of the heavy resonances in the associated process was also discussed. We emphasize that the relevant derivations, prescriptions and strategy are rather general and can be straightforwardly extended to other event topologies, even with invisible particles (see e.g., Ref. [105]).
As a proof of principle, we have applied this general method to study the distinction between the seesaw models based on the Standard Model and the Left-Right symmetric model. Once there is statistically significant evidence for such a jj signal, we can pinpoint the diagram(s) responsible for this signal by using our kinematic endpoint method and also measure the masses of the new particles involved in this process. This can be a powerful way to study the origin of neutrino masses at the LHC and beyond. Along this line, we examined the signal sensitivity for some well-motivated channels at the LHC and a 100 TeV future collider, and found that a W R gauge boson mass up to ∼ 5.5 TeV with 300 fb −1 data at √ s = 14 TeV LHC or up to ∼ 32 TeV with 1 ab −1 data at √ s = 100 TeV collider. This has important consequences for other new physics observables, such as neutrinoless double beta decay, lepton flavor violation and even matter-antimatter asymmetry [128,129].

A Invariant mass distributions
We illustrate the invariant mass distributions for the regions discussed in Section 3 in Figure 12. The relevant distributions for region R 7 are not shown here because all kinematic endpoints (in squared mass) are simply given by either s or s/2, hence less informative than the others. Among the eight possible invariant mass variables (2.2)-(2.4), two are trivial, viz., m jj always peaks at m W R and m jj peaks at m W for some cases; hence they are not considered here. For any event, the decay process of an on-shell resonance is performed through pure phase space. The events are generated with √ s = 14 TeV for all regions. As clearly mentioned in Section 3, the mass of particle C merely determines the overall scale (for the on-shell C-initiated processes), while all the details are completely governed by the mass ratios R ij . In this sense, the exact values of all mass parameters are irrelevant. Every distribution is plotted in the unit of m C which is fixed to be 1 TeV. Figure 12. Invariant mass distributions in regions R 1(1) (first two panels in the 1st row), R 1(2) (last two panels in the 1st row), R 2(1) (first two panels in the 2nd row), R 2(2) (last two panels in the 2nd row), R 3 (first two panels in the 3rd row), R 4 (last two panels in the 3rd row), R 5 (first two panels in the 4th row), and R 6 (last two panels in the 4th row). In each pair of panels, the first one demonstrates the 2-body invariant mass variables such as m (red dashed), m > j (blue dot-dashed), and m < j (green solid) while the second one demonstrates the 3-body invariant mass variables such as m j (red dashed), m > jj (blue dot-dashed), and m < jj (green solid). Black dashed lines denote the relevant theoretical kinematic endpoints formulated in Section 3.