$\Delta L \ge 4$ processes

We discuss the experimental prospects for observing processes which violate lepton number ($\Delta L$) in four units (or more). First, we reconsider neutrinoless quadruple beta decay, deriving a model independent and very conservative lower limit on its half-life of the order of $10^{41}$ ys for $^{150}$Nd. This renders quadruple beta decay unobservable for any feasible experiment. We then turn to a more general discussion of different possible low-energy processes with values of $\Delta L \ge 4$. A simple operator analysis leads to rather pessimistic conclusions about the observability at low-energy experiments in all cases we study. However, the situation looks much brighter for accelerator experiments. For two example models with $\Delta L=4$ and another one with $\Delta L=5$, we show how the LHC or a hypothetical future pp-collider, such as the FCC, could probe multi-lepton number violating operators at the TeV scale.


Introduction
So far, no lepton (L) nor baryon (B) number violating process has been observed experimentally. However, there are good reasons to believe that neither of these quantities are actually conserved. In fact, even within the standard model (SM), non-perturbative effects such as the sphaleron [1] violate both B and L. More phenomenologically, also the observed baryon asymmetry of the universe points to the existence of B violation at some point in its early history.
From the viewpoint of standard model effective theory, one can build non-renormalizable operators which violate B and L [2,3]. The lowest dimensional operator of this kind, the Weinberg operator, appears at dimension 5 (d = 5) and violates L by two units. This operator generates Majorana neutrino mass terms which can be experimentally probed by the observation of neutrinoless double beta decay, 0νββ: (A, Z) → (A, Z ± 2) + 2e ∓ (for recent reviews see for example [4,5]). Next, at d = 6, one finds ∆B = ∆L = 1 operators. These operators cause proton decay in modes such as the famous p → e + π 0 . This and other two-body nucleon decays are well-known to arise in a variety of models, most notably in Grand Unified Theories -see [6,7] and references contained therein.
The gauge structure of the SM and its field content is such that ∆L = 2n+∆B for all non-renormalizable operators (n being an integer). Thus, for example, even ∆L is associated to even ∆B, so no proton decay mode with ∆L = 2 can exist. However, starting at d = 9 one finds ∆L = 3 operators, associated to ∆B = 1. Also, operators relevant for ∆B = 2 processes, such as neutron-antineutron oscillations, appear first at d = 9. One would naively assume that the rates for ∆L ≥ 3 (or also ∆B ≥ 2) processes are necessarily much smaller than those corresponding to the lower dimensional non-renormalizable operators. However, this may not be the case, and in fact it is possible that ∆B, ∆L = 1 processes are forbidden altogether. This is exactly what happens in the Standard Model, since sphalerons are ∆B = ∆L = 3 transitions (thus, sphalerons cannot destroy protons). 1 There is also the possibility that beyond the SM there exists some (so far unknown) symmetry such that lepton number and/or baryon number can only be created or destroyed in larger multiples. For ∆L = 3 this has been recently discussed in [11], see also [12]. In that case, for example, standard proton decay modes are absent and one is left with p → e + νν, π 0 e + νν, e − ννπ + π + and, more interestingly, e + e + e + π − π − . As noted above, these processes are induced by d = 9 or higher operators, which implies that the proton decay rate is suppressed by many powers of the new physics scale Λ. Consequently, one can have Λ ∼ TeV, making it possible for colliders to probe this hypothesis [11].
In this paper, we will discuss ∆L = 4 processes (addressing also the possibility of having ∆L ≥ 5), and analyse whether any of these processes can possibly be observed in the foreseeable future. Specific example models will be constructed where, due to the presence of some symmetry, operators involving less leptons are forbidden.
We start by noting that all ∆L = 4 operators must have dimension greater or equal 10 (see table 1 for examples). Note that at d = 10 there is just one unique ∆L = 4 operator. 2 At d = 12 there are already 8 different operators (including derivatives), so we show only some examples for d ≥ 12. 3 Some operators require more than one generation of fermions, such as the operator at d = 10. Adding two derivatives to this operator, a very similar operator can be realized with just one generation of leptons at d = 12. We will discuss these two operators in more detail in sections 2 and 4. The second example at d = 12 simply exchanges one Higgs from the operator at d = 10 for two quarksQū c , equivalent to a Yukawa interaction. In this way, higher dimensional operators with less Higgs fields can be constructed. This is important, as will be seen in sections 2 and 3, for experiments at low energies. Also at d = 12 we can find an operator involving only fermions (third example in table 1). However, for this particular operator the final state in any process will always involve neutrinos, making it impossible to tag the lepton number experimentally. 1 It is interesting to note that it has long been believed that sphaleron transitions are unobservable experimentally. However, some recent papers on sphaleron rates at accelerators have come to much more optimistic conclusions [8][9][10]. 2 We use the word operator to denote a combination of fields, independently of the number of contractions. In the case of L i L j L k L l HHHH (subscripts denote the lepton flavors), there are two contractions: the square of the Weinberg operator, , and another one O ijkl . Nevertheless, it is not necessary to include in the Lagrangian this last one because it is related to the O contraction with the ijkl indices permuted. Furthermore, note that there are several identities among the entries of the tensor O ijkl , so a total of six couplings/numbers encode all possible L i L j L k L l HHHH interactions (for example, if three or more lepton indices are the same, the operator is identically zero). 3 The complete list up to d = 22, compiled with the Sym2Int program [13], can be found at renatofonseca.net/dl4OpList.php.
We also show in the table two examples of ∆L = 4 operators with d = 15, both of which violate baryon number. The first one has two L fields, thus the SU (2) L contractions again lead to final states involving a neutrino. We will be more interested in the other operator, with four charged leptons. In fact, its realization/decomposition necessarily involves new coloured fields which can be searched at the LHC, as we will point out in section 4.

∆B Operator dimension
Example operators 15ū cūcūcūcūcdcēcēc LL,ū cūcūcūcūcūcēcēcēcēc Table 1: Some examples of ∆L = 4 operators. Note that the lowest dimensional operator at d = 10 is unique. All odd-dimensional operators with d < 23 and ∆L = 4 violate also B, and the smallest ones appear at d = 15.
As far as we know, the only ∆L = 4 process treated in some detail in the literature is neutrinoless quadruple beta decay (0ν4β), having been discussed for the first time in [14]. A simple power counting suggests that the decay rate associated to 0ν4β will be extremely small if all particles which mediate this process are heavy [15]. However, the analysis in [14,15] leaves open the possibility of having an observable 0ν4β decay rate, if several of the mediator particles are very light (here "light" implies masses of the order of the nuclear Fermi scale, i.e. O(0.1) GeV). Thus, in section 2 we calculate a very conservative lower limit on the 0ν4β decay lifetime, based on collider searches for charged bosons. We find that this minimum lifetime is around 20 orders of magnitude larger than the current experimental limit [16], rendering 0ν4β virtually unobservable.
In section 3 we extend the discussion to other ∆L = 4, 5, 6, · · · low-energy processes. In all cases our simple rate estimate is far below experimental sensibilities. The most optimistic case, dinucleon decay of the form (A, Z) → (A − 2, Z − 2) + 2π − + 4e + , is expected to be at least some 8 orders of magnitude beyond the current Super-Kamiokande sensitivity. Thus it seems impossible for low-energy experiments to test lepton number violation in 4 or more units.
However, the prospects of observing ∆L ≥ 4 processes at colliders are good, as we discuss in section 4. Indeed, because operator dimensionality becomes irrelevant at energies comparable to the new physics scale, the LHC is in a good position to probe the creation or destruction of groups of 4 or more charged leptons. We construct two example models for ∆L = 4 (and one for ∆L = 5) and calculate production cross sections for the different particles in these models at pp-colliders (the LHC and a hypothetical √ s = 100 TeV collider). From these cross sections we derive lepton number violating event rates and estimate the scales up to which pp-colliders can test such kind of models. We then end the paper with a summary.

Minimum lifetime of neutrinoless quadruple beta decay
The possibility of observing ∆L = 4 processes via neutrinoless quadruple beta decay (0ν4β) was put forward in [14]. Three candidate nuclei for the decay (A, Z) → (A, Z + 4) + 4e − were identified: by far the most promising one is 150 Nd, which could decay into 150 Gd plus four electrons with a total kinetic energy of Q = 2.08 MeV. 4 With 36.6g of the isotope, the NEMO-3 detector is well suited to measure this decay provided that it happens at a reasonable rate. Recently, based on an exposure of 0.19 kg·y of 150 Nd, the collaboration reported a lower limit of (1.1 − 3.2) × 10 21 y on the half-life for this particular process, at a 90% confidence level [16]. The range in this number is explained mostly by the fact that, so far, no reliable theoretical calculation of the single electron spectra has been done for this process. Nevertheless, as can be seen, the result is fairly insensitive to such details.
Given the non-observation of 0ν2β decay events so far, one has to wonder if 0ν4β decays will ever be observed. If the only contribution to the latter process is the conversion of a pair of neutrons into protons and two electrons, twice repeated (see fig. 1), then one would expect the approximate relation between the lifetimes of double and quadruple beta decay without neutrinos. 5 In this expression q ∼ 100 MeV stands for the typical momentum transfer in a nucleus. Using the experimental lower limit on τ 0ν2β for various nuclei of the order of 10 26 years [18,19], we can extract the lower bound τ 0ν4β 10 88 years. However, this estimate assumes that the main contribution to neutrinoless quadruple beta decay comes from two virtual double beta decays. This does not need to be the case; indeed, as previously mentioned, it is possible to forbid entirely ∆L = 2 processes and still have those with ∆L = 4.
In the following we will argue that an important constraint on neutrinoless quadruple beta decay can be derived using data from collider searches of charged bosons. However, before proceeding, let us consider first what is the expected value of the 0ν4β decay rate from simple power counting (see also [14,15]). We start by noting that there are 12 fermions involved, so after electroweak symmetry breaking the relevant operator has dimension 18: Here, κ is just some unspecified dimensionless coefficient. On the other hand, in the final state there are four electrons and a nucleus which stays essentially at rest, so the 0ν4β decay rate will depend on the (3 × 4 − 1)th power of the available kinetic energy Q. With these considerations, inserting some numerical factors for the multi-body kinematics, we obtain the formula Clearly, this is a very rough estimate for the life-time of the process, with an uncertainty of a few orders of magnitude. Nevertheless, given the largeness of the numbers involved, it is good enough for the following discussion.
To obtain the lowest possible τ 0ν4β , ideally one would need large, i.e. order O(1), couplings (κ ∼ 1) and light mediator masses (Λ i TeV). Note that for very light particles, their mass becomes irrelevant when compared to the typical momentum transfer in the nucleus. Hence the best case scenario is Λ min i ∼ q. One can see easily that, in the limit where all Λ i are of the order of 1 GeV or lower, it becomes conceivable to have τ 0ν4β 150 Nd smaller than 10 26 years. However, we will now show that the 0ν4β decay diagram contains at least 4 propagators of charged bosons (scalar or vector). The lightest particle, that can play this role, is the W boson, hence the most optimistic scenario achievable is This corresponds to roughly one 0ν4β decay per year for a mass of ∼ 10 17 kg of neodymium, hence the observation of this process would be extremely challenging even in the most optimistic scenario. Independent on any concrete model, four or more of charged bosons propagators are needed for the following reason. Conservation of fermion number implies that the 12 external fermions which make up the 0ν4β operator must be arranged in six currents J i . Out of the six possible pairings (uu, ud, ue, dd, de and ee) none is electrically neutral. Thus all six J i currents must exchange charge through scalar or vector bosons. (It is also conceivable that this charge exchange occurs through some intermediary vertex, yet this scenario does not lead to a minimal number of internal propagators.) It is straightforward to see that the most economical setup is the one, where pairs of currents J i with opposite electric charged are connected by a single charged boson, in which case only 3 such propagators are needed. However, the only fermion billinears with opposite charges are dd and de, and clearly one does not obtain the 0ν4β operator with three copies of these fermions. Hence, a minimum of 4 charged bosons are needed and the lower limit in eq. (3) applies.
The six possible currents need to be coupled to scalars or vectors. These are charged or doubly charged particles, leptoquarks or diquarks. If one were to construct a loop model for 0ν4β decay and introduce also some new exotic fermions, more exotic scalars/vectors could, in principle, also appear. All of these states, however, necessarily couple to standard model fermions and thus can be searched at accelerators such as LEP and LHC. Again, considering the large numbers involved in eqs. (3) and (4), a very rough argument suffices for our purpose. Thus, we only quote that LEP data rules out any electrically charged boson, decaying to SM fermions, below roughly 100 GeV [20], The lower bound in eq. (4) on the neutrinoless quadruple beta decay lifetime is therefore unavoidable. However, it is also very conservative. If one were to construct lower bounds individually for the 6 possible currents, (much) larger limits could be derived for the different individual cases. Instead, we will now try to see, whether it is actually possible to approach this bound, by considering some particularly promising models.
We have established already that 0ν4β decay is necessarily suppressed by the heavy mass of at least 4 charged bosons and, in order not to further reduce the decay rate, one should avoid colored particles in internal lines. 6 It is easy to check that, out of the 135 distinct ways of partitioning the 0ν4β operator in fermion billinears/currents J i , only (ud) (ud) (ud) (ud) (ee) (ee) makes it possible to have all internal bosons colorless. It is equally simple to arrive at the conclusion that, for this particular partition of the 0ν4β operator, the number of neutral and colorless internal bosons will be minimal (just 4) if and only if the corresponding diagram can be split into two halves, each with fermion billinears (ud) (ud) (ee), connected by a neutral boson. A particularly interesting way of building such diagram is the one shown in fig. (2), where only SM fields plus a neutral scalar are used [14]. In the central part of the diagram, it is clear that 4 neutrinos are created from nothing, hence the neutrino-scalar interaction is critical for the violation of lepton number. Under the full Standard Model symmetry, and setting aside a complication which we shall mention later, this ∆L = 4 operator must be of the form LLLLHHHH (or perhaps higher dimensional). Crucially, this effective interaction does not need to be suppressed by the smallness of neutrino mass; indeed, the operator LLHH might even by absent. One way of generating the LLLLHHHH interaction is by adding to the Standard Model an extra SU (2) L scalar triplet ∆ with one unit of hypercharge, as well as a real scalar singlet σ with no gauge interactions. The most general Lagrangian with these two new fields violates lepton number in two units, given the simultaneous existence of the couplings LL∆ and HH∆ * . 7 It is possible to construct a model which violates lepton number in units of 4 only, and to do so we retain just the following interactions: It is straightforward to arrange scalar potential parameters such that only the Higgs doublet acquires a non-zero vacuum expectation value, in which case there is no spontaneous breaking of the Z 4 symmetry. As such, the Weinberg operator is not generated, but LLLLHHHH is (see fig. (3)). Assuming without loss of generality that the coupling κ is real, σ mixes only with Re ∆ 0 , resulting in two mass eigenstates with masses m 1,2 and a mixing angle θ which can be controlled at will. In particular, it is possible to (a) make one of these states very light, (b) have maximal mixing sin θ, and at the same time (c) keep all other new scalars (Im ∆ 0 , ∆ + and ∆ ++ ) arbitrarily heavy. Neutrino masses can be accounted for by introducing right-handed neutrinos, or by the breaking of this Z 4 symmetry, for example by introducing a very small HH∆ * coupling. Likewise, the 0ν2β decay rate can be zero (absence of HH∆ * ) or very small (small ∆L = 2 coupling).  The complication mentioned earlier is that there is no 1 Λ 2 ν e,L ν e,L ν e,L ν e,L local operator: given that the spinor ν e,L has two components only, the Pauli exclusion principle forbids point interactions with more than two ν e,L . Instead, the four-neutrino interaction must be of the form where Λ is a mass associated to the particle(s) mediating the effective interaction. We used the fact that the distance between the two interaction points x and y is inversely proportional to Λ. In practice, this means that there is an extra suppression factor (momentum/Λ) 2 due to Pauli blocking on top of the expected 1/Λ 2 suppression for this 4-fermion interaction 8 , and that the relevant operator is ∂∂LLLLHHHH instead of LLLLHHHH.
As pointed out in [14], the uncertainty in the invisible decay width of the Z boson as measured indirectly at LEP [22][23][24][25] constrains the mass of the scalar boson in fig. (2) to be heavier than roughly Λ ∼ 20 GeV. In this latter process, the momentum p is of the order of the Z mass, hence p ∼ Λ and there is no Pauli blocking. However, the blocking is present in the 0ν4β decay diagram with 4 W 's, since p ∼ 100 MeV is much smaller than the mass Λ of the neutral scalar. As a consequence, the 0ν4β decay lifetime is even larger than what has been assumed previously [14,15]; a simple order of magnitude calculation yields τ 4W 0ν4β ∼ 10 69 years, which is very far from the limit given in eq. (4). A possible way of avoiding this particular suppression factor is by introducing right-handed neutrinos ν R and W R gauge bosons such that lepton number is violated by a 4-fermion interaction 1 Λ 2 ν e,L ν e,L ν e,R ν e,R . This effective operator can be generated with the scalars ∆ and σ mentioned earlier, keeping in mind that there is now a ν e,R ν e,R σ interaction. Nevertheless, given the multi-TeV LHC mass limits on the new gauge bosons [26,27], the two diagrams imply a similar 0ν4β decay life-time: Finally, we will make some comments about the possible realizations of the 0ν4β decay operator in eq. (2) with Standard Model fields (see also [28]). Due to Pauli's exclusion principle, fermion fields evaluated at the same space-time point anti-commute, so for a generic 4-spinor we can only have operators of the form Ψ n [· · · ] with n ≤ 4 (if there are no derivatives). 9 Furthermore, if we expand the Dirac indices of such an operator with n = 4, the only non-zero term must be proportional This means that the Pauli exclusion principle severely restricts local 0ν4β operators to the unique form and higher dimensional 0ν (2m) β decay operators, with m > 2, are forbidden entirely unless they have derivatives. Note, that quarks have 3 colors, hence a similar issue arises for operators with more than 6 quarks of the same charge and chirality. 8 This kind of amplitude suppression can equally appear in operators with no fermions. For example, consider a global SU (2) symmetry under which the scalars φ = (φ 1 , φ 2 ) T and φ = φ 1 , φ 2 T transform as doublets, and S is a scalar singlet. The local operator φ (x) φ (x) S (x) S (x) is identically zero, but φ (x) φ (y) S (x) S (y) is not, hence it is possible to collide two S's and create a φ 1 plus a φ 2 . This is done via two diagrams in which φ 1 and φ 2 are exchanged in the t and u channels. In the limit where the mass Λ of φ 1,2 is much larger the momentum transferred p, the total amplitude varies with 1/Λ 2 × (p/Λ) 2 , and not 1/Λ 2 . 9 If only one chirality X = R, L is involved, then Ψ n X [· · · ] with n ≤ 2 are the only possibilities. For example, an electronneutrino mass term ν e,L ν e,L is allowed, but interactions of the form ν e,L ν e,L ν e,L [· · · ] or ν e,L ν e,L ν e,L ν e,L are not.
In this sense, quadruple beta decay (with or without neutrinos) is a borderline case between allowed and excluded local processes. An interesting consequence is that there are only three 0ν4β operators of minimal dimension (=18): This should be contrasted with the 0ν2β decay operators of dimension 9, of which there are six. On the other hand, the diagram shown above with 4 W 's ( fig. (2)) and the equivalent one with 2W 's + 2W R 's, have dimensions 24 and 20 respectively:

Other low energy processes with ∆L ≥ 4 involving charged leptons
We now move on to a brief discussion of other lepton number violating processes, involving low-energies, with four or more charged leptons. Rough estimates for their rates will be given in the following. We stress that for an experimental proof of L violation, final states should not contain neutrinos.
In table 2 we give the lowest dimension at which a given ∆L = 0 operator can appear, together with some examples. (Note that, since we are interested here in low energy processes, neither Higgs nor gauge bosons can appear as final states.) The table starts with the ∆L = 1 operators which induce the standard proton decay modes (hence ∆B = 1), and these are followed by the ∆L = 2 operator associated to neutrinoless double beta decay. With larger ∆L, the dimension of the operators keeps rising and at some point one expects that the rates of associated low-energy processes becomes too small to be observed. We now discuss briefly this point, by focusing on the most promising processes. 24 e 6 u 8 d 2 · · · · · · · · · It is important to distinguish those cases where baryon number is violated from those scenarios where where ∆B = 0. This is simply due to the fact that the available energy in ∆B = 0 processes is fixed by the nucleon mass, of the order of ∼GeV, while kinetic energy of the charged leptons is much smaller (∼MeV) in the ∆B = 0 case.
Let us consider first the latter case, ∆B = 0. This implies immediately that ∆L must be an even number. The relevant processes are then 0ν (2n) β with n > 2. We will discuss only β − decays, since for quadruple beta decays it has been shown already in [14] that the positron emission or electron capture processes are even more hopeless, due to their smaller Q-values. For 0ν (2n) β with n > 2 the same observation applies.
The process 0ν6β is induced by an operator with 18 fermions, hence the decay width Γ is suppressed at least by a factor Q 17 q 29 /Λ 46 relative to the nucleon Fermi momentum q ∼ 100 MeV. Moreover, we note that at most 4 electrons can be at a single point x, therefore the operators for these 0ν (2n) β decays require at least two derivatives, and consequently the decay width is suppressed by 4 more powers of q/Λ, compared to the simple-minded estimate quoted above. It is also straightforward to check that at least 6 electrically charged bosons are needed to mediate the process hence, the same logic as discussed for neutrinoless quadruple beta decay applies.
Finally, kinematically 0ν6β and larger is only allowed for neutron-rich nuclides which are far from the valley of stability, hence these isotopes will have a very short half-life. The longest-lived isotope seems to be Let us now turn to processes where baryons are destroyed and hence, the available energy is much larger. Here, we consider the cases (∆L, ∆B) = (4, ±2) and (∆L, ∆B) = (5, ±1); violation of lepton or baryon number in greater quantities will be associated with even bigger minimum life-times.
The lowest dimensional operator with 4 charged leptons (after electroweak symmetry breaking) is eeeeuuuuuu/Λ 11 (see table 2). It leads, for example, to diproton decay: (A, Z) → (A − 2, Z − 2) + 2π − + 4e + . In the most optimistic scenario, this operator can be built in such a way that only 7 powers of Λ correspond to the mass of mediators with colour (Λ C ), while the remaining 4 powers of Λ are related to the mass of fields with electroweak interactions only (Λ EW ). Hence, using the values Λ EW ≈ 200 GeV and Λ C ≈ 2 TeV. The prefactor 10 −13 takes care of the fact that this is a 6-body decay. Super-Kamiokande has searched for other diproton decay modes, imposing limits of the order of 10 32 years on the associated life-times [29]. It seems therefore very hard to probe ∆L = 4 processes at low energies, even for those cases where baryon number is violated. Note, however, that due to the much larger energy release, the gap between the experimental sensitivity and the most optimistic expectation is only 8 orders of magnitude, compared to the (minimum of) 20 orders found for quadruple beta decay. For ∆L = 5 (or larger) the decay rates are necessarily even more suppressed. For ∆L = 5, the lowest dimensional operators have d = 21, as shown in table 2, therefore the decay widths will be suppressed by 34 powers of ∼ m p /Λ when compared to m p . In summary, we conclude that observation of charged lepton number violation in four or more units in low energy experiments will be impossible in the foreseeable future.

Probing lepton number violation at colliders
We now turn to a discussion of probing models with multiple lepton number violation at colliders. We do not aim at a full, systematic analysis of all possibilities. Instead, we will discuss two simple models with ∆L = 4 and then present one example for a model with ∆L = 5. Models which lead to different ∆L = 4 operators at low energy or models with larger ∆L violation can be easily constructed following the same principles that we use in our examples.
Our first model is inspired by the discussion on neutrinoless quadruple beta decay in section 2. In this model, called model-I below, we add only two new fields to the standard model. Both are scalars: (i) ∆ = S 1,3,1,−2 and (ii) T = S 1,3,0,−2 . Here, the subscripts stand for the transformation properties or charge under the SM gauge group and lepton number, SU (3) C × SU (2) L × U (1) Y , L (also we use an S for scalars and, later on, an F for 2-component Weyl spinors). The only change with respect to the model discussed in section 2 is that we have replaced the singlet σ = S 1,1,0,−2 with the Y = 0 SU (2) L triplet field T .
As before, we enforce a Z 4 (L) symmetry which ensures that leptons can only be created or destroyed in groups of 4. The Lagrangian of the model is: 10 Note that this Lagrangian is also U (1) B invariant. In other words, baryon number is preserved, hence processes such as dinucleon decay are completely absent in this model. Note that ∆ is the same field that appears in the type-II seesaw mechanism. However, our symmetry forbids the term H∆H, which in seesaw type-II is the source of ∆L = 2. This implies that for m 2 ∆ ≥ 0, in our model there is no induced vev for ∆ 0 and thus no Majorana neutrino mass term.
We have not written generation indices in eq. (16). In general Y ∆ is a complex symmetric (3,3) matrix. All terms in the Lagrangian, with the exception of those proportional to T 2 or T 4 , conserve lepton number. For the phenomenology discussed below it is important that m 2 T violates ∆L in 4 units. The term proportional to λ HH∆T leads to mixing between the neutral and singly charged components in T and ∆ after electro-weak symmetry breaking (EWSB), as well as to a mass splitting between the CP-even and CP-odd components in ∆ 0 . Thus, after EWSB the model has two new neutral CP-even scalars S 0 1,2 and two singly charged scalars S ± 1,2 plus one neutral CP-odd scalar, A 0 , and one doubly charged scalar, ∆ ±± . In our numerical calculations we always diagonalize all mass matrices and consider mass eigenstates correctly. However, in the following discussion we simply use T 0 and ∆ 0 for S 0 2 and S 0 1 , respectively. (And similarly for the singly charged states T ± and ∆ ± .) This is done only for the clarity of the discussion; it does not affect any of our conclusions. Note that for typical choices of masses m 2 T and m 2 ∆ above (500 GeV) 2 mixing between the different states will be small unless m 2 T m 2 ∆ . For the numerical calculations shown below, we have implemented the model into SARAH [30,31]. The implementation is then used to generate SPheno code [32,33] for the numerical generation of spectra. The UFO model files generated by SARAH are used for cross section and decay calculations with MadGraph [34]. Both, ∆ and T can be produced with sizeable rates at colliders. Fig. (4) shows cross sections for the most important production modes in pp-colliders for two values of √ s: To the left √ s = 13 TeV, to the 10 We use the notation ∆ ≡ ∆ ++ ∆ + TeV. The numerically largest cross section is found for ∆ ±± pair production. However, as we will discuss below, for the observation of ∆L = 4 processes the interesting production modes are pp → (T 0 + T + ) and pp → (T 0 + T − ). Both processes proceed through an off-shell W -diagram, see fig. (5). The cross section for pp → (T 0 + T + ) is larger than for pp → (T 0 + T − ), reflecting the fact that the initial state is positively charged.
We first discuss the decays of ∆. The different components of ∆ decay according to ∆ ±± → l ± α l ± β , ∆ ± → ν α l ± β and ∆ 0 → ν α ν β with 100% branching ratio, when summed over α and β. Since cross sections are largest and background lowest for ∆ ±± , the most stringent constraints on m ∆ come from searches for ∆ ±± . Both ATLAS [35] and CMS [36] have searched for doubly charged scalars decaying to charged leptons. Limits depend quite strongly on the flavour of the charged leptons. CMS [36] gives limits as low as m ∆ ±± 535 GeV for a ∆ ++ decaying with 100 % to pairs of taus, while limits are in the range of (800 − 820) GeV, if the ∆ ++ decays only to electrons or muons. ATLAS [35], on the other hand, has established lower limits on m ∆ ±± of roughly (600 − 800) GeV, for branching ratios to either electrons or muons in the range of (0.2 − 1). We will therefore use two choices of m ∆ in our numerical examples below, namely, m ∆ = 0.6 TeV and m ∆ = 1 TeV. The former is allowed only for a ∆ with coupling mostly to τ 's, while the latter is currently unconstrained. Note, however, that with the predicted L = 3/ab for the high luminosity LHC (HL-LHC) m ∆ in excess of m ∆ = 1 TeV will be probed, while for a √ s = 100 TeV collider we estimate from fig. (4) that up to m ∆ ∼ 5 TeV could be tested with L = 3/ab, in agreement with the numbers quoted in [37].
The coupling λ HH∆T was chosen λ HH∆T = 0.1 in this example, while the entries in Y ∆ where arbitrarily put to be smaller than O(0.1). With these choices, decays to purely leptonic final states are negligible and therefore not shown.
For the observation of a ∆L = 4 processes, we need to produce the hyperchargeless triplet T . We therefore now turn to a discussion of the decays of T 0 and T ± . First of all, note that all decay rates for these particles will be proportional to the coupling λ HH∆T . This is due to the fact that this term is the only one linear in T allowed by Z(4) L . Mixing between T 0 and T ± with ∆ 0 and ∆ ± induce two-body decays for these states into leptonic final states. However, these will always involve neutrinos and thus are not useful to establish experimentally lepton number violation. More important are then decays of T 0 and T ± to ∆ ±± and gauge bosons. Fig. (5) shows the most important Feynman diagrams. Apart from T + → W − + ∆ ++ , T + can decay to W + + A 0 , W + + ∆ 0 , 11 as well as ∆ + + h and ∆ + + Z 0 . The branching ratio for T + → W − + ∆ ++ is always close to 25%, if the mixing between T + and ∆ + is small. Similarly, T 0 can decay to a number of final states involving gauge (and Higgs) bosons. Fig. (6) shows the most important branching ratios for T 0 decays as a function of m T . As the plot shows, for small values of m T the two-body decays T 0 → A 0 + Z 0 (and T 0 → ∆ 0 + h) have the largest branching ratios. However, for large values of m T , the decay mode T 0 → ∆ ±± + 2W ∓ becomes dominant. T 0 decays to both ∆ ++ and ∆ −− with equal rates. That both decays have the same rate can be understood as a mass insertion in the decay of T 0 , see fig. (5). Recall that m 2 T violates L by 4 units. At a pp collider one can therefore have lepton number violating (conserving) final states from the chain With the cross sections, see fig. (4), and the decay branching ratios we can then calculate the number of ∆L = 4 events as a function of m T (and m ∆ ). We have performed this exercise for the LHC and find that, taking into account the lower limit on m ∆ , there will be less than 1 event even after 3/ab of data has been taken. However, the prospects look much brighter at an hypothetical √ s = 100 TeV collider, see fig. (7). In this figure we show the maximal number of ∆L = 4 events attainable in model-I for two choices of m ∆ and two different values for the expected luminosity. For L = 30/ab more than 10 ∆L = 4 events could be found up to roughly m T 6 TeV. Thus, for our first model we conclude that while the LHC can extend the search for the different particles in model-I to above 1 TeV, observation of ∆L = 4 seems not to be possible at the LHC. At a √ s = 100 TeV collider ∆L = 4 could be discovered in this model up to a scale of roughly 6 TeV.
The negative conclusions for the LHC can be simply understood from the fact that model-I contains no new coloured states. As shown in table (2), the smallest dimensional operator generating four charged leptons is 15-dimensional: e 4 u 6 . We therefore choose to implement it in our second example, model-II. Any model leading to this operator will necessarily involve beyond-the-SM coloured fields.
Model-II introduces three new states, two fermions O = F 8,1,0,−2 and D c = F3 ,1,1/3,1 (together with its vector-partnerD c = F 3,1,−1/3,3 ) and one scalar, S d = S 3,1,−1/3,1 . By enforcing once again Z 4 (L) invariance, we obtain this time a Lagrangian with an enlarged accidental symmetry group G SM × U (1) 2B−L . It might not be immediately obvious that this latter group contains Z 4 (L), but this is nevertheless true. Indeed, the Lorentz and Standard Model group G SM force all operators with SM fields to be Z 2 (B − L) invariant, hence we may write L = 2n + B for some integer n. Together with L = 2B, it is then quite easy to see that L and B are forced to be multiples of 4 and 2, respectively. Crucially, unlike in model-I, due to the U (1) 2B−L symmetry it is not possible to break lepton number without breaking baryon number as well. As such, one can have small/unobservable dinucleon decays rates, but neutrinoless quadruple beta decay is strictly forbidden.
The Lagrangian contains the following terms: We have calculated the pair production cross sections for the new particles in our model-II using again MadGraph [34]. The results are shown in fig. (4). Again, the plot to the left is for the LHC, the one on the right is calculated for √ s = 100 TeV. In all cases one expects that gluon-gluon fusion gives the largest contribution to the cross section, see fig. (9). The largest cross section is found for the fermionic octet. More than 10 events from O-pair production are expected in L = 3/ab for octet masses up to m O 3 TeV. A √ s = 100 TeV collider would be able to collect more than 10 events for octet masses up to m O 15.5 (18.5) TeV for L = 3/ab (30/ab).  In model-II, since O is an electrically neutral state, it will decay with equal branching ratios to D c 1/3 + S d,−1/3 and D c −1/3 + S d,1/3 . This implies that ∆L = 4 final states 4l + + 6j will have the same rate as the ∆L = 0 final states 2l + + 2l − + 6j. Thus, the production cross section (event number) of O-pair production gives directly the limit on the scale up to which ∆L = 4 can be tested in model-II. It is straight-forward to use the ideas discussed above to construct also models, which lead necessarily to larger ∆L. We will discuss only one example with ∆L = 5. This model-III introduces five new states.
We will need two copies of S d , to which we assign different lepton numbers, S L=1 d = S 3,1,−1/3,1 and S L=0 d = S 3,1,−1/3,0 . The vector-like down quarks also now come in two copies. We have D c L=2 = F3 ,1,1/3,2 , its vector partner isD c L=3 = F 3,1,−1/3,3 , and also D c L=1 = F3 ,1,1/3,1 with its vector partnerD c L=4 = F 3,1,−1/3,4 . Finally, the model also has the fermionic vector octet, O = F 8,1,0,2 andŌ = F 8,1,0,3 . With this lepton number charges, we then enforce a Z 5 (L) symmetry. Just as with model-II, there is a bigger, accidental symmetry group in this model, U (1) 5B−3L . In other words, for each group of 5 leptons created, 3 new baryons should appear as well, and for this reason the proton is completely stable in this model-III. ∆L = 5 processes at the LHC can then occur through diagrams such as the example shown in fig. (10), whereD c L=3 is pair produced via gluon fusion. Note that the decay chains of bothD c L=3 and D c L=3 * end with the same number of SM fermions: Seven. One can assign the source of lepton number violation in this diagram to the mass termD c L=3 D c L=2 . If all other particles in the diagram are lighter than D c L=2 and all couplings the same order, ∆L = 5 and ∆L = 0 final states from these decay chains will have similar rates. 12 Nevertheless, note that even if theD c L=3 D c L=2 mass term was switched off, lepton number would still be broken; in fact, B and L conservation would be restored only if the vector masses of D c L=1 and O were absent as well.
We can estimate the mass reach of the LHC to test this kind of diagram from the cross sections shown in fig. (8). We estimate that more than 10 events (before cuts) would remain for masses of D c L=2 below 2.7 TeV in L = 3/ab. At a √ s = 100 TeV collider more than 10 events would occur for m D c below 13.3 (15.5) TeV in L = 3/ab (30/ab). Thus, there is amble parameter space that could be tested in future colliders even for models with ∆L > 4.. Many different models of this kind can be readily constructed.

Summary
Given the current experimental situation, the total number of leptons L might be a conserved quantity. Standard probes for ∆L = 0 are proton decay (∆L = 1, 3) and neutrinoless double beta decay (∆L = 2) experiments. However, neither have found any signal so far. It is therefore possible that L is violated only in larger multiplicities, i.e. three, four or more units.
In this context, we have discussed that the decay of a nucleus into four electrons and no neutrinos will likely never be observed. This and other ∆L = 4 low energy processes where the leptons are electrically charged must necessarily be mediated by several heavy bosons, hence their amplitudes are severely suppressed, given the current accelerator constraints. We have calculated a very conservative lower limit on the half-life for neutrinoless quadruple beta decay, which is 20 orders of magnitude larger than current experimental sensitivities. The same conclusion is valid if lepton number is broken in five or more units. The least pessimistic scenario that we have found is the one where baryon number is also violated. In this case two protons in a nucleus could decay into 4 positrons plus pions, for example. The energy scale in this process is set by the proton mass, but even so, the high dimensionality of the operators involved implies that proton decay experiments would need to increase their exposure by at least 8 orders of magnitude before meaningful constraints could be derived experimentally.
Colliders, on the other hand, can explore the possibility that lepton number is violated in four or more units. The reason behind this observation is rather simple: Even though these processes involve many exotic particles, if the energy of the collider exceeds the mass of those exotics, the suppression associated with the high dimensionality of these ∆L ≥ 4 operators disappears. In this work, we presented two models for ∆L = 4 and one for ∆L = 5, which make use of this idea. We have calculated cross sections for the LHC and a possible future √ s = 100 TeV collider, estimating the rates for ∆L = 4 (and ∆L = 5) processes. Naturally, for hadron colliders the expectations are highest for models which contain coloured particles. In this case the LHC (the √ s = 100 TeV collider) could probe ∆L = 4 up to scales of roughly 3 TeV (18 TeV). We expect that the high multiplicity of these events associated with ∆L ≥ 4 should make them virtually background free, giving a rather spectacular signal.