Discriminating the HTM and MLRSM models in collider studies via doubly charged Higgs boson pair production and the subsequent leptonic decays

We present a case study for the doubly charged Higgs bosons $H^{\pm\pm}$ pair production in $e^+e^-$ and $pp$ colliders with their subsequent decays to four charged leptons. We consider the Higgs Triplet Model (HTM) not restricted by the custodial symmetry and the Minimal Left-Right Symmetric Model (MLRSM). The models include scalar triplets with different complexity of scalar potentials and, due to experimental restrictions, completely different scales of non-standard triplet vacuum expectation values. In both models, a doubly charged Higgs boson $H^{\pm\pm}$ can acquire a mass of hundreds of gigaelectronvolts, which can be probed at HL-LHC, future $e^+e^-$, and hadron colliders. We take into account a comprehensive set of constraints on the parameters of both models coming from neutrino oscillations, LHC, $e^+e^-$ and low-energy lepton flavour violating data and assume the same mass of $H^{\pm\pm}$. Our finding is that the $H^{\pm\pm}$ pair production in lepton and hadron colliders is comparable in both models, though more pronounced in MLRSM. We show that the decay branching ratios can be different within both models, leading to distinguishable four lepton signals and that the strongest are $4\mu$ events yielded by MLRSM. Typically we find that MLRSM signals are one order of magnitude larger that in HTM. For example, the $pp \to 4\mu$ MLRSM signal for 1 TeV $H^{\pm \pm}$ mass results in a clearly detectable significance of $S \simeq 11$ for HL-LHC and $S \simeq 290$ for FCC-hh. Finally we provide quantitative predictions for the dilepton invariant mass distributions and lepton separations which help to identify non-standard signals.


Discovery of the first Higgs scalar boson at LHC and the Standard Model
The spectacular discovery of the chargeless Higgs particle (H 0 ) at the LHC [1][2][3] is consistent with the prediction of the Standard Model (SM), confirming the basic concept of the spontaneous symmetry breaking mechanism and elementary particle mass generation. Observed H 0 decays into gauge boson particles W + W − and ZZ [4][5][6] fits beautifully into this picture. Similarly, determination of ttH 0 couplings in gluon fusion [7][8][9] and ttH 0 production [10] confirm Higgs boson role in fermion mass generation. With gathered statistics, we know more and more about this particle, namely its decay rate to γγ [11,12], spin-parity which is dominantly J P = 0 + [4,5,[13][14][15]. Also mass suppressed decay rate to muon pairs when comparing to top pairs is evident [12,16]. Yet another spectacular success of the LHC physics is a clear discovery that the Higgs boson decays to the third generation of fermions, namely to the pairs of τ leptons and b-quarks. Especially determination of the Yukawa Higgs boson coupling to b-quarks is tricky, as though this channel amounts to about 60% of Higgs boson decays, the QCD b-quark background is overwhelming [17]. The story of the Higgs boson studies continues. Very recently the measurements of the Higgs boson's properties have reached a new stage in precision through detection of a rare decay mode where the Higgs boson decays into two muons [18,19]. Aiming at sub-percent precision for Higgs boson decays, quantitative tests of the SM for Higgs boson couplings need further scrutinization in studies at HL-LHC and future Higgs factories. These also include the investigation of the Higgs boson self-coupling [20].

Searching for new scalar bosons at future colliders and a choice of tested models
Detection of the Standard Model scalar particle does not preclude validity of more elaborated physical scenarios with extended scalar sectors. The simplest extensions beyond the SM doublet scalar multiplet include their copies, like the two Higgs doublet model [21], supersymmetric extensions of the SM [22,23] or, stepping up in this construction, scenarios with triplet scalar representations either in their supersymmetric [24,25] or non-supersymmetric versions [26][27][28]. Here we will consider the latter. There are many possibilities for triplet representations, depending on the hypercharge Y ≡ 2(Q − T 3 ) [29][30][31][32]. We will explore the simplest one which involves doubly charged Higgs fields in the triplet representation with hypercharge Y = 2, the Higgs Triplet Model (HTM) [33]. For that, we will not assume any special symmetries or constructions [34,35], so that v ∆ , the triplet vacuum expectation value (VEV) will be extremely tiny, at the scale of electronvolts which makes experiments more challenging. We will also consider a much more complex model where the Standard Model SU (2) × U (1) gauge symmetry is extended by an additional SU (2) group, the so-called minimal left-right symmetric model (MLRSM) [28,[36][37][38][39][40]. Thus we consider a setting where both the HTM and the MLRSM models include doubly charged Higgs bosons.
HTM received a considerable amount of attention recently [41][42][43][44][45][46][47][48][49][50][51][52][53][54]. This model when confronted with experimental data, features a strong restriction in which v ∆ is very small, O (1) (GeV), or below. Here, in particular, we concentrate on the cases where v ∆ is of the order of neutrino masses. Then the triplet Yukawa couplings will be of the O (1) order and H ±± decays dominantly into the same-sign dilepton channel. In this case, the LHC direct search bound on the doubly charged scalar mass, m H ±± 850 GeV [55] applies. At the same time, the constraints from different lepton flavor violating (LFV) processes and non-universality of leptonic couplings start to weigh in. There is thus a direct relationship among the triplet VEV v ∆ , neutrino masses, their mixing and doubly charged Higgs couplings. That is why the production and decays of H ±± scalars at high energies depend substantially on the oscillation data and limits on LFV processes in HTM.
On the other hand, in MLRSM the dominating non-standard effects in phenomenological studies are connected with the right-handed breaking scale v R which affects the couplings and masses of a wide set of non-standard heavy particles of spin 0,1,1/2 present in the model. Low-energy precision  = 700 GeV. The scalar potential parameters, fields and relations for masses are defined in the Appendix, Eq. 7.1 and Eq. 7.18. We identify h and H 0 0 as the SM Higgs boson (H 0 ).
As we can see in Tab. 1, there are more scalar fields in MLRSM than in HTM. Then any detectable signal connected with a neutral, singly or doubly charged Higgs bosons which are present in MLRSM but are not present in HTM would be in favour of MLRSM. However, we should note that though MLRSM is very rich in particle content and non-standard interactions, despite enormous theoretical and experimental efforts over last several decades, what we get so far are exclusion limits on the parameters of this model. All experimental data considered so far gives no indication for neutral, singly or doubly charged scalars, extra neutral heavy leptons, extra gauge bosons. So, the starting point is actually the same: we do not know if and which BSM model is realized in nature and we are still looking for a first experimental indication towards any non-standard signals in one or another model. As the models' parameters are already severely constrained, we have to consider very rare processes and hence faint signals. We focus on the cleanest BSM colliders signal connected with doubly charged scalars: their pair production and subsequent decays (correlations between the same-sign leptons in the final state originated from H ±± decays). To leave no stone unturned, we will focus especially on the case in which two doubly charged Higgs bosons in MLRSM have the same masses as otherwise the second scalar H ±± 2 connected with right-handed triplet in MLRSM (see Appendix for fields definitions) would help to discriminate between both models in favor of MLRSM. A case with different H ±± 1 and H ±± 2 masses in e + e − CLIC center of mass energies will be shortly discussed in Section 5.5, non-degenerate mass cases for hadron colliders have been discussed already in [88]. For the same masses of H ±± 1 and H ±± 2 , the production rates are higher in MLRSM than in HTM. We will see what is the contribution of H ±± 2 against H ±± 1 in production processes at lepton and hadron colliders, in case of the same doubly charged boson masses and how these contributions change with center of mass energy. As we will see, there are scenarios with model parameters where the difference of signals in both models can be further enhanced by studying leptonic branching ratios of doubly charged Higgs bosons and kinematic cuts.
Fixing the scalar mass spectrum lets us take a first look at the production processes e + e − → H ++ (1,2) H −− (1,2) and pp → H ++ (1,2) H −− (1,2) for HTM and MLRSM. This will bring us to the discussion of the importance of the neutrino sector. Fig. 1 and Fig. 2 show classes of Feynman diagrams for the H ±± pair production in e + e − collisions in both models. There are s-channel diagrams mediated by neutral gauge bosons Z and γ and Higgs bosons, and the t-channel diagram. Due to experimental restrictions discussed in the next two sections, the contributions coming from the s-channel diagrams are comparable to the off-resonance regions, and the resonance regions for the considered center of mass energies and masses lie away from the allowed region of parameters (see Figs. 7,8). It gives possibility to discuss how the t-channel diagrams in Fig. 1 and Fig. 2 affects the process.
As schematically depicted in the figures the relevant H ±± − l ∓ − l ∓ vertices come from Yukawa couplings. e + e − In HTM the Yukawa term (7.5) with neutrino fields generates Majorana masses This term also contains the H ±± − l ∓ − l ∓ vertex which leads to lepton flavor violation. We can diagonalize the neutrino mass matrix by U as follows [89] The matrix U relates the mass eigenstates |ν i through a superposition of the flavor states |ν : |ν i = U T |ν , so it is directly connected with the PMNS matrix (3.7) and the exact relation between them is U * = V PMNS . Now we can write the Yukawa couplings as a function of the PMNS matrix and the masses of neutrinos. From Eq. (2.2), Y can be written in the following form We discuss the parametrization of V PMNS and the employed range of the oscillation parameters in Section 3.1.
The Y coupling depends on v ∆ , neutrino masses and oscillation parameters. From perturbativity, Y 2 ≤ 4π. Apart from this restriction, there are stringent limits on Y coming from various experimental data discussed in the next section.
In MLRSM the t-channel with the H ±± −l i −l j vertex is inversely proportional to v R . We assume, see section 7.3, vanishing off-diagonal couplings. In this case the vertex is As we can see, the coupling Eq. (2.3) in HTM can be enhanced in the case of small values of v ∆ → 0. However, it is at the same time proportional to the light neutrino masses. The analogous coupling H ±± −l i −l j in MLRSM is related to the heavy neutrino masses and v R , which are limited by, e.g. bounds on heavy gauge boson masses, see section 4. In the next two sections, we will consider details of the considered models to find the allowed space of the models' physical parameters, including the neutrino sector, which as seen enters the considered processes with very different light and heavy masses.
In general, it would be tempting to find a way to show when the processes of doubly charged Higgs boson pair production decouple from the neutrino masses. Though such relations are a feature of considered models, if the signals which we predict in both models would not fit experimental data, this would be a sign for another mechanism that takes place. For MLRSM and HTM the basic neutrino mass mechanisms are the seesaw type-I and type-II, respectively.
The H ±± pair can be produced in the proton-proton collider via photon, Z boson and neutral scalar particles in the s-channel, see Fig. 3. As will be discussed in section 5.1, due to existing experimental constraints, also here the production process is very similar in both models. What will bring the difference are doubly charged Higgs boson decays which lead to the four charged lepton final signals. To discuss it properly, in the next two sections we will present relevant experimental constraints on the models' parameters.

The HTM model and the relevant experimental parameter constraints
The Yukawa term (7.5) includes the H ± − l − ν and considered in the last section H ±± − l − l vertices. They can contribute to several LFV processes like radiative decay of charged leptons l i → l j γ, three body decay of charged leptons l i → l j l k l l , µ-to-e conversion in nuclei µN → eN * . We show the contributing diagrams for HTM in Fig. 4(a)-(e). In Fig. 4(f) we include the muoniumantimuonium conversion µ − e + → µ + e − . Corresponding limits on the H ±± parameters are gathered in Tab. 2. Tab. 3 gathers relevant SM processes: Møller scattering e − e − → e − e − , Bhabha scattering e + e − → e + e − , the anomalous magnetic moment of the muon (g − 2) µ from which also useful limits on the H ±± parameters are derived. These processes have been discussed in the context of HTM in many works [33,[90][91][92][93][94][95][96]). The branching ratios (BR) in Table 2 depend on the charged scalar masses and the Yukawa couplings Y , defining the allowed space of mass and coupling parameters for charged scalars. The radiative decays l i → l j γ and the µ-to-e conversion process mediated by doubly and singly charged scalar bosons originate from the following Lagrangian [97] Branching ratios depend on the form factors A L and A R , which actual form depends on Higgs scalar contributions to the considered processes. For the doubly charged scalar there are four relevant diagrams as shown in Fig. 4 (a) and (b). The amplitude for H ±± for the first two diagrams Fig. 4 (a), at the leading order of the doubly charged scalar mass is µ is a mass parameter introduced in dimensional regularization, = 4 − D and D is dimension. The contribution from other two diagrams Fig. 4(b) mediated by the doubly charged scalar boson is By adding (3.2) and (3.3) we can see that the final contribution is finite and after doing some algebra the contribution of the doubly charged scalar form factors can be written in a compact way In a similar way one can compute the contributions from diagrams mediated by singly charged scalar bosons and the total amplitude in HTM can be written as (m e P L iσ µν q ν + m µ P R iσ µν q ν ). (3.5) Matching Eq. (3.5) to Eq. (3.1), we can extract the form of the A L (A R ) form factors and compute the analytic formulas for the radiative decays and µ-to-e conversion processes. The final analytic formula for considered LFV processes are gathered in the Appendix. If the massive neutrinos couple to leptons and are of Majorana type, the lepton number can be violated by two units, ∆L = 2. This leads to the neutrinoless double beta decay ββ0ν process [98,99] and as it has not been observed so far, it puts a constraint on the model parameters. This process has been analyzed within HTM in [100], the non-standard contribution is negligibly small.
Above we have discussed LFV processes which have not been observed so far, leading to stringent bounds on BSM physics and parameters. Useful information about limits on the BSM physics can be also obtained by exploring observed SM processes and analysing experimental results and SM predictions. One of suchfinite processes is the Bhabha scattering present in electron-positron collisions. It serves as a calibration method for colliders as it is a QED dominated t-channel process, see Section II and Figs. 1-2 in [101]. The LEP data [102] sets a lower limit on H ±± , namely Y ≥ 10 −7 (to ensure H ±± decay before entering the detector).
Another SM experiment which seems to provide a promising signature of the BSM physics is the observed excess in the anomalous magnetic moment of the muon (g − 2) µ . There is a lasting discrepancy of more than 3σ in the measurement of (g − 2) µ with the corresponding SM value [103]. At present the deviation, as given by PDG, is [104]: The experimental limits for Bhabha, Møller and (g − 2) µ SM processes are collected in Tab. 3. Charged scalars can contribute to the (g − 2) µ at the one-loop level. There are many studies of BSM contribution to (g − 2) µ in the literature. The contribution from a doubly charged Higgs boson to (g − 2) µ is discussed in [105] and in the context of HTM in [91]. The diagrams mediated by singly and doubly charged scalars contributions to (g − 2) µ are given by Fig. 4 (a) and (b) where both l i , l j are µ (muons). Contributions of singly and doubly charged scalar bosons to (g − 2) µ amounts a negative number [57] and (g − 2) µ anomaly is hard to explain by H ±± . However, it is worth mentioning that the observed anomaly is an open problem as still some discrepancies among different low-energy experiments exist [103].    3) we can see that the H ±± − l − l couplings depend on the neutrino oscillation parameters, neutrinos hierarchy and the lightest neutrino mass. Details of studies for the HTM model are thus very sensitive to the neutrino oscillation data, as discussed already in [33] and [54,86,120]. In our analysis the following, standard parametrization of the V PMNS matrix is used: where s ij and c ij denotes sin(θ ij ) and cos(θ ij ), respectively. Tab. 4 shows global neutrino fits at the 2σ C.L. for neutrino parameters which are used in present analysis for two mass orderings, defined as: Normal mass hierarchy: Inverted mass hierarchy:  Table 4: Neutrino oscillation data, notations as in [121]. ∆m 2 ij = m 2 i − m 2 j . Depending on the hierarchy, for atmospheric nutrino oscillations either ∆m 2 3l = ∆m 2 31 > 0 (NH) or ∆m 2 3l = ∆m 2 32 < 0 (IH).

Process
Concerning the Dirac CP-phase δ CP , the global fits indicate preference for its non-zero values. Recent T2K results confirm this tendency and considered by us 2σ range of the Dirac phase covers well the best fit values given in [122]. In analysis we are chosing δ CP data as given in Tab.4. There is no direct limit on the Majorana phases α 1 , α 2 . However, in some studies there are predictions using the neutrinoless double beta decay, e.g., see [123]. There is no bound on the individual masses of neutrinos from the oscillation data. Therefore, the lightest neutrino mass m ν0 is a free parameter and other two masses are determined through (3.8). Also, there are limits on the sum of three neutrino masses from different experiments: from the tritium decay [124] or neutrinoless double beta decay [125], the sharpest limit comes from astrophysics and cosmology [126] Σ ≡ 3 i=i m νi ≤ 0.23 eV. (3.9) These limits set the upper bound on the lightest neutrino mass [127,128], present experimental data gives m ν0 = 0.071 eV, NH, 0.066 eV, IH. (3.10)

The triplet VEV v ∆ and the ρ-parameter
As we mentioned in the introduction, the additional scalar triplet contributes to the ρ parameter. It can be defined either through a relation among massive SM gauge bosons Z and W and Weinberg mixing angle or relations among gauge couplings [129]. In HTM, at the tree level, ρ can be written as [130]: (3.11) The experimental limit on the ρ parameter [131]: ρ exp = 1.00037 ± 0.00023, (3.12) put the upper bound on the triplet VEV v ∆ .
Taking [132,133] where G F is Fermi coupling constant 1.1663787(6) × 10 −5 GeV −2 [134], we get v ∆ ≤ 1.7 GeV, (3.13) for ρ exp , within 2σ deviations. Let us note that the limit on v ∆ can not be obtained for ρ exp within 2σ deviation. It is connected with the fact that relation (3.12) has sense only for ρ exp ≤ 1, otherwise v ∆ comes out to be a complex number. We will see in the following section that other low experimental data are more important, lowering down the scale of v ∆ in an unambiguous way to the (sub)electronvolt level.

Relation between v ∆ and doubly charged scalar particles parameters in the light of low and high energy experimental limits
In this section, we analyze bounds on the triplet VEV v ∆ from low and high energy experiments discussed earlier (see Tables 2 and 3). Fig. 5 shows excluded regions in the plane of v ∆ and M H ±± parameters' space based on current limits on branching ratios (for both NH and IH scenarios) for various LFV processes and (g − 2) µ . The analytic formulas for the relevant quantities are collected in the Appendix. In analysis we consider 2σ range of neutrino oscillation parameters, Tab. 4, Majorana phases α 1 and α 2 are varied in the full range (0,2π). We vary lightest neutrino mass m ν0 keeping the Σ (sum of neutrino masses) limit (3.9) for both inverted and normal hierarchies. We assume degenerate mass for charged scalar bosons, M H ±± = M H ± , and vary them from ∼ 500 GeV to 1000 GeV (M H ±± 470 is already excluded by the LHC, see section 5 and a discussion around Tab. 9). The shaded regions in Fig. 5 are excluded from LFV and muon (g − 2)µ limits.
(a) (b) Figure 5: Plots for v ∆ vs m H ±± using normal and inverted hierarchy data. Shaded regions correspond to the exclusion limits coming from LFV bounds for current data and future sensitivity expectations. The neutrino oscillation data are taken in the 2σ range. In general, precision of future experiments (see Tab. 2) will allow to get one order of magnitude better limits on v ∆ .
We use different colors to show the exclusion from individual LFV processes: radiative decay of leptons (green), three body decay of leptons(red), µ-to-e conversion (blue), and (g − 2) µ (violet).
We can see that the most stringent limit is due to three body decays l i → l j l k l k specifically the µ → eee process. We do not find any significant difference between two neutrino mass hierarchy scenarios, but for low neutrino masses the radiative decay µ → eγ starts to play an important role in normal hierarchy case (see Tab. 5). Bounds coming from scattering processes or muonium to antimuonium conversion are at least one order of magnitude smaller than those obtain through (g − 2) µ calculation and are not included in the above plots.
In Table 5, we collect the lower limits of v ∆ in eV for different values of Majorana phases and lightest neutrino mass assuming m H ±± = 700 GeV. The process which put the strongest limit is written below the numerical values. For our further analysis and the HTM benchmark point we take v ∆ = 15 eV.

The MLRSM model and relevant experimental constraints on its parameters
We consider a left-right symmetric model based on the SU [28,[36][37][38]40] in its most restricted form, so-called Minimal Left-Right Symmetric Model (MLRSM) which contains a bidoublet Φ and two (left and right) triplets ∆ L,R [29,38,39,88] see the Appendix for details.

Constraints on MLRSM model parameters and the triplet VEV v R
The heavy sector of the model is triggered by VEV v R connected with the Higgs triplet ∆ R . All new gauge and scalar bosons are proportional to v R , and v R κ, where κ is a VEV related to the scale of the SM spontaneous symmetry breaking and to the SM gauge bosons W 1 , Z 1 , κ 246 GeV, see Eq. (7.17).
Using the relation between the heavy charged gauge boson mass and the SU we can find the parameter space for v R and heavy neutrino masses. In the last few years the LHC has constrained the possible v R scale very much by exploiting different channels where W 2 plays a crucial role, e.g., W 2 decays to two jets [135], two jets and two leptons [136] and top-bottom quarks [137]. Altogether, the following bounds on M W2 have been obtained: (i) ATLAS -3.  [136], assuming that SU (2) R gauge coupling g L equals the SU (2) L coupling g R . These bounds can be relaxed without such an assumption [139][140][141][142]. The CMS experimental data based on the pp → lljj process are presented as the M W2 − M N exclusion plots, see Fig. 6 in [136] and Fig. 7 in [143]. For convenience we repeat them here in Fig. 6. We use these data, and analogous data from the ATLAS collaboration [144], leading to restrictions on the t-channel in Fig. 8 and final signals presented in Section 5.5. a) b) c) As assumed in Fig. 6 we take M W2 ≥ M N and a correlation between the masses which are proportional to v R [56,145]. By passing, let us note that most of the experimental LHC analyses are based on simplified scenarios where heavy neutrinos are mass degenerate with diagonal mixings and where CP-violating effects are not taken into account. However, CP-parities of neutrinos can be different leading to destructive interference effects, relaxing limits on the v R scale, see [146][147][148].
A simultaneous fit to the SM low energy charged and neutral currents set a rather weak bound on M W2 , so v R , namely M W2 > 715 GeV [104,146]. However, there are additional restrictions for model's parameters coming from radiative corrections. As far as one loop corrections are concerned and additional precision constraints on MLRSM parameters, there are very few studies based on LR models, i.e., [39,129,146,149] (MLRSM model), the other papers are: [150,151] (limits on W 2 mass coming from the K L − K S process (finite box diagrams, renormalization not required)), [152] (LEP physics), [153][154][155][156][157] (process b → sγ). Some interesting results are discussed also in papers [158,159] where the problem of decoupling of heavy scalar particles in low energy processes has been discussed. In [160] it has been shown that low-energy radiative corrections shrink non-standard parameters to very small regions, due to correlations among gauge bosons, scalars and heavy neutrino massses, though still there is a freedom connected among others with unknown scale v R . We assume v R and windows of possible masses of heavy MLRSM particles allowed by low energy analysis [160][161][162][163][164].
Apart from experimental limits, due to tree-unitarity and flavor changing neutral currents (FCNC) constraints, the scalar potential parameter α 3 in Eq. (7.18) is restricted and masses of neutral Higgs bosons H 0 0 , A 0 should be greater than 10 TeV. The lowest limit on v R scale is 1.3÷6.5 TeV [56], depending on the mass scale of FCNC Higgs bosons [165]. Such a relatively low (TeV) scale of the heavy sector is theoretically possible, even if gauge unification (GUT) is demanded, for a discussion, see [166] and [167].

5.1
The H ±± pair production at e + e − and pp colliders As discussed in previous sections, we assume M H ±± (1,2) = 700 GeV. This value will be further justified when the H ±± decay branching ratios are discussed in next sections. Therefore, for substantial H ±± pair production in e + e − collisions, we need the centre mass energy √ s above 1 TeV. As discussed in Introduction such energies for e + e − colliders are planned presently only at CLIC. Numerical results for √ s = 1.5 TeV are gathered in Fig. 7 and Fig. 8 for HTM and MLRSM, respectively.
vertex is proportional to the left-handed triplet VEV v L which is set to zero to preserve the ρ-parameter [169].
As discussed in Section 2, the t-channel in HTM contains the e − l − H ±± vertex inversely proportional to v ∆ in Eq. (2.3), this diagram becomes dominant for small v ∆ . However, it appears that the region where the t-channel can dominate is ruled out by the low energy data and Tab. 5. The allowed t-channel cross section for e + e − → H ++ H −− is a few orders of magnitude lower than the s-channel, which is equal to 2.4 fb, see a solid horizontal line in Fig. 7. As it is shown, regardless of the choice of the neutrino parameters, the whole region where the t-channel is not negligible is excluded.
cross section in MLRSM, see Fig. 2 depends on the righthanded triplet VEV v R and heavy neutrino masses. The allowed space for M W2 − M N parameters has been considered in section 4.1 and is based on limits on the heavy neutrino masses taken from the LHC CMS and ATLAS data for the pp → lljj process [136,143,144]. This process is a collider analogue of the neutrinoless double beta decay mediated by a heavy charged boson, heavy Majorana neutrinos, and cross-sections depend strongly on masses and CP-parities of heavy neutrinos [147]. As we have in disposal CMS and ATLAS results, in calculations we assume M W2 > M N with the same CP-parities of heavy neutrinos. In Fig. 8 we vary the M W2 mass from 600 GeV to 5.5 TeV and the heavy neutrino mass up to 4.8 TeV and take the best fit expected values for the LHC exclusion data.
The production through the t-channel is constrained by the Yukawa coupling, Eq. Y ee (2.4). We assume perturbativity of the coupling with h M 1 we get the relation between v R and heavy neutrino masses. Since the LHC exclusion plots assume M N < M W2 , this condition is fulfilled automatically for the considered parameter space.
The most strict limits comes from the Bhabha and Møller processes, see Fig. 9, the doubly charged scalar particles can contribute there. In Tab. 6 we gathered region of physical masses for heavy neutrinos which arise from the discussed low energy LFV constraints.
(a) (b) Figure 9: The e + e − → e + e − (Bhabha) and e − e − → e − e − (Møller) processes at the lowest order with doubly charged Higgs bosons. As we can see in Fig. 8, the t-channel gray parts of the plotted lines above the long-dashed "Bhabha, Møller" line assigned with cross × and plus + symbols might dominate within the whole region of the v R parameter tested by LHC. However, adding the discussed Yukawa constraints on H ±± couplings gathered in Tab. 2 and Tab. 3, this region is eliminated (corresponding allowed t-channel contributions with red and blue parts of the plotted CMS and ATLAS lines are thickened in Fig. 8). As the Bhabha and Møller processes constraint the t-channel contribution to be below 0.3 fb, altogether with the LHC constraints, it results in a much smaller contribution than the s-channel contribution and the interference effect is small: The total cross section σ tot practically corresponds with the s-channel. Even though the mass M Z2 is a function of v R (M Z2 0.78 v R ), the higher resonances are suppressed since the small center mass energy is too small to observe them. For larger v R values, we are outside the s-resonance for √ s = 1.5 TeV and s-channel contributions are flat and small. For instance, for v R = 6(15) TeV which will be used as reference values in next sections for four lepton final state analysis, and which correspond to M Z2 = 4.7(11.7) TeV, σ s 4.6 fb. The limits from the muon (g − 2) and the µ + e − → µ − e + process are also taken into account, since the corresponding diagrams contain the f ee and f µµ couplings, but they play no significant role. The (g − 2) µ process restricts the f µµ coupling, see the Appendix. It affects heavy neutrino mass bounds and for further calculations we assume that maximum M N2 = 5 TeV, what is safe for considered values of v R (6 and 15 TeV). Unlike in the HTM case, the LFV processes do not restrict further the results because we assume the LFV vertices to be negligible with no light-heavy neutrino mixings (see Section 7.3). Taking into account the above constraints, the maximal cross section at the t-channel is σ t ∼ 0.3 fb.
All non-standard heavy particle masses are related to the vacuum expectation value of the righthanded triplet, see Appendix 7.2 and Eqs. (7.21)-(7.29). As discussed in [56], the combined effects of relevant Higgs potential parameters and Higgs bosons responsible for FCNC limits regulate the lower limits of heavy gauge boson masses. In Fig. 8 we put only low-energy limits on v R coming from (g − 2) µ and FCNC. We indicate v R ∼ 3.5 TeV, which by considering the Higgs boson mass spectrum Eqs. (7.21)-(7.27) is the minimal v R for FCNC Higgs masses of A 0 1 , H 0 1 scalars at the level of O(10) TeV, and the minimal allowed M H 0 3 for α 3 scalar parameter to be less than 16. The mass limit for A 0 1 , H 0 1 at the level of 10 TeV is the lowest limit on FCNC Higgs boson masses [165], one of the strongest limits has been obtained in [154] (M A 0 1 ,H 0 1 ≥ 50 TeV). We can see that there are various estimates of the v R scale, see also Fig. 6. Apart from the dijet LHC strong limits, there are searches in the one jet and one lepton signal category [138,170] as well as off-shell W 2 and Z 2 channels [171,172]. All these studies confirm that it is not natural to expect v R to scale below 3.5 TeV. For these reasons, as pp studies at HL-LHC or future FCC-hh or CEPC colliders offer investigation of heavy BSM states at higher scales, in next sections for pp phenomenological studies we assume v R scale and MLRSM mass benchmarks corresponding to higher v R values at the level of 6 and 15 TeV.
In Tab Process  Let us proceed to the hadron colliders and pair production of H ±± Higgs bosons. Basic treelevel diagrams for considered models are given in Fig. 3.   14 TeV 100 TeV pair production at HL-LHC c.m. energies while for FCC-hh/CEPC option, Z 2 -channel starts to be important. Due to the shown differences between MLRSM and HTM models we can expect a higher number of events for 4-lepton final states in MLRSM when masses of doubly charged Higgs bosons are the same. However, it does not have to be the case as the final results depend strongly on branching ratios which we will consider in the next section.
The QCD contributions to the doubly charged Higgs boson pair production increase the cross section at the NLO level. The role of the QCD effects in the hadronic processes of H ±± pair production has been considered in [54]. A similar situation with positive contribution of QCD at the NLO and higher levels has been observed also for other processes in models which include triplet Higgs bosons and heavy neutral leptons [54,[173][174][175][176]. The corresponding k-factors (which measure ratios of higher order QCD effects to the tree level cross section) do not change considerably with the H ±± mass and centre of mass energies, k-factor ∈ (1.15 ÷ 1.20). Due to different ratios of H ++ 1 and H ++ 2 pair production processes (see Eqs.(5.2)-(5.4) and Tab. 8), for m H ±± = 1 TeV, the k-factor in HTM is 1.15 and is smaller than k-factors in MLRSM, which are 1.6 (1.85) for HL-LHC (FCC-hh/CEPC) centre of mass energies, respectively. There are various QCD contributions at the NLO level to the considered process, in which the s-channels γ/Z 1 /Z 2 dominate over the gluon and photon fusion mechanisms, both for HL-LHC and FCC-hh/CEPC. Concerning potential contributions beyond NLO, the N 3 LL terms are found to be about three times larger than NLO terms. However, this is connected mainly with gluon fusion which is subdominant for the considered H ++ 1 masses in the s-channel [54]. As the doubly charged pair production signals are dominated by the exchange of the SM particles in e + e − collisions (see Tab. 7), differences between doubly charged pair production signals in both models are small. A better estimation of QCD corrections, evaluating the NNLO terms, would resolve expected signals better. In the pp collision case, the production difference between the models for the considered benchmark points is much larger. NLO QCD corrections seem to be enough to discriminate the models, though we should note that the production difference between both models will decrease above v R = 15 TeV, which is the upper limit for the v R value considered in the present work. In scenarios with v R > 15 TeV a knowledge of NNLO QCD corrections will be also useful in the pp collisions. Anticipating final four-lepton results, the above conclusions do not change for the considered benchmark points and kinematic cuts. Namely, ratios of MLRSM (v R = 15 TeV) to HTM four-lepton signals can be as large as 1.7 (34 and 43) at e + e − and pp (HL-LHC and CEPC/FCC-hh) colliders, respectively (see 4µ signals, Tab. 15 and Tab. 16). Then NLO QCD k-factors should be enough to distinguish the HTM and MLRSM signals in pp collisions, unless the v R scale is too large and the Z 2 gauge boson contribution becomes at the NNLO QCD level.
To summarize, the QCD contributions to the considered production processes at the NLO level are substantial in both models and must be taken into account in the analysis. To discriminate both models, evalauation of higher order QCD terms may be needed for higher v R scales.

HTM, a choice of benchmark parameters and H ±± decay scenarios
In HTM the doubly charged scalar has nine possible decay channels, depending on the scalar boson mass (i) H ±± → l i l j , i, j = e, µ, τ , In this paper we focus on the first channel (i) and present a case study for pair production of a doubly charged scalar boson and its subsequent leptonic decays, considered also in [177]. It is a very clean channel which provides a unique signature for colliders signal with a pair of the same sign leptons [88]. Scenarios (iii) and (iv) require non-degenerate masses for charged scalar particles: In Fig. 11 we show a variety of branching ratios as a function of v ∆ for various H ±± decay channels.
On the left plot we show the following decay modes: leptonic (red), W ± gauge bosons (green) and H ± W ± (blue). On right we give a variation of leptonic and pair of gauge boson decay branching ratios for a degenerate mass of H ±± . There are two cases there: the solid line is for M H ±± = M H ± = 700 GeV and the dashed line is for a charged scalar boson mass of 300 GeV (this mass is already excluded by LHC, we left it for comparison with previous work, see Fig. 4 in [178]). The shaded region is connected with the lightest neutrino mass and mass hierarchy, within 2σ  oscillation parameter range. This region does not change the result substantially. We can see that the cross-cut point is shifting with charged scalar boson mass, but in the interesting mass region, the lepton channel dominates till v ∆ reaches values in range of 10 4 ÷ 10 5 eV. In Fig. 11 (a) we take a mass gap M H ±± − M H ± = 80 GeV, in Fig. 11 (b) there is no mass gap and both H ±± → H ± W ± and H ±± → H ± H ± channels are suppressed. It has been shown in [86] and [87] that there are limits on the mass gap |M H ±± − M H ± | in order to preserve the oblique T-parameter, unitarity and potential stability condition. For recent work on vacuum stability conditions of Higgs potentials in various variants of HTM models, see [179]. From electroweak precision data and limits from the h → γγ process [83,85] the dominant contributions are in the degenerate mass case. Therefore only leptonic and W gauge boson decay channels are possible. However, the H ±± − W ∓ − W ∓ vertex is proportional to the triplet VEV v ∆ while the Yukawa coupling in the H ±± − l ∓ − l ∓ vertex is proportional to 1 v∆ , so the lepton channels dominate strongly over the scenario (ii) for the triplet VEV v ∆ < 10 5 eV.
For VEV v ∆ in a range of eV, the cumulative leptonic channel dominates in that region regardless of the neutrino masses and oscillation parameters as well as doubly charged scalar boson masses. So, our final conclusion is that when H ± W ± and H ± H ± channels are suppressed, the leptonic decays dominate for low v ∆ .
The sharpest limit from ATLAS on M H ±± is that the H ±± mass should be larger than 870 GeV for the left-handed triplet doubly charged scalar boson field, assuming the 100% branching ratio for the H ±± → l ± l ± decay (l ± = e ± , µ ± ). However, it is possible to lower down the limit to 450 GeV for a 10% leptonic decay branching ratio (see Fig. 13 d in [55]). On the other hand, the decays into a τ lepton are not considered in the above analysis. In Tab. 9 we present branching ratios for those channels and the result for the ee, eµ and µµ decays, within the ±2σ range of the oscillation parameter space. For other channels including the τ we refer to [180]. The strength of lepton decay channels depends strongly on the neutrino masses, their hierarchies and oscillation parameters. It is possible to find the parameter space where the branching ratio for the particular lepton channel  Table 9: Lowest limits on a mass of the doubly charged scalar boson M H ±± for different branching ratios [55]. We removed data which corresponds to the branching ratio region beyond what has been obtained in Fig. 12 within 2σ range of the neutrino oscillation parameters.
is small regardless v ∆ even if the cumulative lepton channel dominates over the W boson channel (the relative lepton decay contributions Γ(H ±± → l l )/ Γ(H ±± → l i l j ) do not depend on the triplet VEV v ∆ ).
We combine the data both from the LHC limits [55] and neutrino parameters within the ±2σ range given in Tab. 4 and compute the lowest limit on the doubly charged scalar boson mass 2 . In Tab. 9 we removed the BR values that are forbidden due to the neutrino oscillation parameters. Another interesting conclusion from this table is that within the HTM the doubly charged scalar boson cannot be lighter than 473 GeV for the normal neutrino mass scenario (and 518 GeV for the inverted mass hierarchy), see Fig. 12 (a). Finally, the lowest mass limit on M H ±± within HTM is 473.7 GeV for NH and 645.4 GeV for IH with BR(H ±± → ll ) = (Γ(H ±± → e ± e ± + e ± µ ± + µ ± µ ± ))/Γ(H ±± → i,j l ± i l ± j ) ≥ 0.1 and 0.4, respectively, where l i,j = e, µ, τ . The most severe limit at 734 GeV comes from the same sign muon channel when BR is 50%.
In conclusion, when assuming the complete scenarios with H ±± decays to all the leptons, still M H ±± can be relatively light.
Our main aim is to analyse the final four lepton (4l) signals which can be potentially seen at the colliders. The dominant signatures are e + e + e − e − and µ + µ + µ − µ − final states within both HTM and MLRSM models. In MLRSM they are not bounded by the neutrino oscillation parameters since the H ±± 1,2 − l − l vertex is related to the heavy right-handed neutrino masses and parameters, as discussed in section 7.3. Within the HTM model these 4l contributions are restricted by the light neutrino oscillation data. Using branching ratios shown in Fig. 12 we compute two parameter sets (for normal and inverted hierarchy) for which the branching ratio for e ± e ± and µ ± µ ± are the highest. We collect the chosen parameters in Tab. 10. We choose two benchmark masses for the collider analyses: M H ±± = 700 GeV (which can be probed at very high energies in e + e − collision, when available, see section 5.1) and M H ±± = 1000 GeV (this higher mass range can be probed without problems at the HL-LHC and FCC-hh, see Fig. 10). For the e ± e ± decay channel we chose   Fig. 12 (b) and (d) with (M H ±± = 700 GeV) and • (M H ±± = 1000 GeV).

MLRSM, a choice of benchmark parameters and H ±±
1,2 decay scenarios Contributing vertices to the non-leptonic decay channels stem from the kinetic term and scalar potential (see Eqs. 19 and 25 in [39]). Relevant decay modes of doubly charged scalar bosons and respective strength of couplings are gathered in Tab. 11. The emboldened processes in the table dominate for v L = ρ 4 = 0 and ξ → 0 [56,169,181], see the Appendix for details. Apart from the values of vertices, we need to take into account the mass spectrum. To suppress the FCNC processes some of neutral scalar particles have to be heavier than 10 TeV. As a consequence, the mass of H ± 2 should above 10 TeV, (7.21) and (7.27). Therefore we can neglect the H ±± 2 decay to the H ± 2 scalar boson for CLIC and LHC energies. From (7.21) it is easy to find that the triplet VEV should fulfil an inequality: v R > √ 2 10 3 / √ α 3 [GeV]. Because α 3 is a quartic coupling (four-scalar interaction) it contributes to the 2 → 2 scattering and the unitarity condition requires α 3 < 8π [56]. The triplet VEV v R has to be higher than ∼ 2800 GeV that translate to M W2 > 1325 GeV. So we can neglect the doubly charged scalar bosons pair production with the subsequent decay to the heavy gauge boson pair H ±± 2 → W ± 2 + W ± 2 for energies lower than 2M W2 . In Tab. 11 we present the other possible decay channels of H ±± 1,2 and corresponding vertices. Most of them are negligible due to model's consistency [56,169], only the bold decay channels can be sunstantial. The H ±± 1,2 decay to H ± 2 is not possible for CLIC and LHC energies because of the FCNC limits (7.27). Vertices contributing to the H ±± 1,2 decays to W 1 , W 2 can be large and are included in analysis leading to final four lepton signals.
Regarding H ±± 1 , its decay to H ± 1 + W ± 1 is limited by Higgs potential parameters and, as proved analytically in [182], the allowed split ∆M can not exceed value 65.3 GeV. We choose the benchmark points for v R = 6 TeV and 15 TeV. The first value falls in energy range of LHC with pp → W 2 → lN l → llW * 2 → llqq , assuming that M Ni < M W2 . Corresponding experimental results can be found in [136,143]. We assume that the doubly charged scalar masses are degenerate and choose two benchmark points: 700 GeV and 1000 GeV.  Table 11: Doubly charged scalar boson decay channels to scalar and gauge bosons in MLRSM. We have listed all possible vertices, thickening the dominating processes assuming that the left triplet VEV v L is equal to zero and keeping in mind experimental limits on the W 1 − W 2 mixing angle ξ < 10 −2 [134,146] and setting the ρ 4 parameter to zero [56,169]. The leptonic decays are analysed separately.
100% also in a case of H ±± 1 . Here the situation is different than in HTM where upper bounds for H ±± branching ratios are given. As discussed in Sections 2 and 3, neutrino Yukawa couplings can be rewritten in terms of oscillation parameters and v ∆ and experimental data restricts possible branching ratios in a substantial way. In addition, depending on the branching ratios, the lowest limit on the mass of a doubly charged scalar can be obtained. However, in the context of MLRSM, leptonic branching ratios for H ±± depend in addition on v R scale and heavy neutrino masses and couplings. This freedom makes it possible to reach full leptonic decays for M H ±± , as given in Tab. 12 and Tab. 13.  are due to right-handed leptonic couplings as analysed in [55].
Tab. 13 shows chosen, allowed values of heavy neutrino masses for given maximal branching ratios. They are consistent with assumption M W2 ≥ M N discussed in section 4, and a correlation between the masses which are proportional to v R [56,145]. In section 5.2 we have obtained the lowest limits for M H ±± as a function of the lightest neutrino mass for a given H ±± branching ratio, arguing that M H ±± at the level of 700 GeV is still possible within HTM. The lowest limits on masses of H ±± 1,2 Higgs bosons have been obtained in [169] by analyzing restrictions on the scalar potential. Before we present final results we will discuss the SM background for the considered leptonic final states.   [136,143]. The heavy neutrino masses for v R = 6 TeV fulfill the low energy constraints given in Tab. 6. M N1 is mostly restricted by the Møller scattering, while M N2 is bounded by (g − 2) µ .

The four leptons background in pp and e + e − collisions
We are interested at estimation of the Standard Model background for pp → l + l − l + l − and e + e − → l + l − l + l − processes, where l ± = e ± , µ ± . The four leptons production at LHC is discussed in [88,183]. The most relevant processes which contribute to the background are tt(Z/γ * ) and (Z/γ * )(Z/γ * ) production. To optimize the collider non-standard effects (decreasing SM tri-and four-lepton SM background and reducing the efficiency of misidentification of b-jets as leptons), we use the following criteria and selection cuts C1. Lepton identification criteria: transverse momentum p T ≥ 10 GeV, pseudorapidity |η| < 2.5.
C6. Hadronic activity cut -within the cone of radius 0.2 around the lepton the hadronic activity should fulfill the inequality: The results are gathered in Tab. 14. For the e + e − collision we consider scattering and annihilation channels with photon radiation, (Z/γ * )(Z/γ * ) production and multiperipheral processes in Fig. 13. The most relevant are diagrams b) and d). For √ s = 1500 GeV we get σ = 4.465 fb before and σ = 0.415 fb after applying the cuts defined above.   [186], accordingly, for tt(Z/γ * ) it is k = 1.6 [187], for (Z/γ * )(Z/γ * ) it is k = 1.5 [188]. Cross section values are given in fb.

Final four lepton signals within the HTM and MLRSM models, a comparison
The final signals depends on subsequent H ±± decays (→ 4e, 4µ, 2e2µ) and suitable kinematic cuts. In the HTM model we take benchmark points for the model connected with maximal 4e and 4µ signals as given in Fig. 12    The 4-lepton signals obtained for the e + e − case are gathered in Tab. 15. In section 5.4 we defined the kinematic cuts which maximise the 4-lepton signals. With assumed total luminosity, we can see that the SM background is comfortable small for muons and the maximal 4µ signal's prediction in HTM can be significant, which is not true in the case of electrons. The difference is enhanced by assumed detector efficiency for electrons (muons), see the cut C2 in section 5.4. For MLRSM chosen parameters in Tab. 13 and v R = 6 TeV the signals are small when compared to the SM background, especially for electrons. For muons a signal is ∼ 3 times smaller than in HTM. However, for v R = 15 TeV the signals for muons detection can be larger in MLRSM, since for that value of v R , M N and M W2 values of parameters lie outside the region examined by the LHC (Fig  6). In this case independent, maximal branching ratios for H ±± 1,2 → e ± e ± and H ±± 1,2 → µ ± µ ± can reach 100%, (Tab. 13), which is not possible for HTM (Tab.10).
It should be noted that the e + e − → 4l results in MLRSM depends strongly on interference effects and the chosen heavy neutrino parameters as the LHC exclusion data affects directly the t-channel contributions. In fact, comparing the HTM results with the MLRSM results for v R = 6 TeV, we can see that the 4l signals can be larger in HTM where the t-channel is negligible for all allowed parameters space (Fig. 7), while in the MLRSM model the t-channel effects can still be large and comparable to the s-channel contributions (Fig. 8). However, in both models the signals are much below the SM background level.
The maximal significance value S ≡ S / √ S + B where S and B are the total number of signal and background events is S = 14 for 4µ signals in MLRSM with v R =15 TeV. For HTM, S = 11 for both NH and IH neutrino mass scenarios and the 4µ signal. The goal for HL-LHC is to deliver about L = 0.25 ab −1 = 250 f b −1 per year with the aim of integrating a total luminosity in the range of 3 to 4.5 ab −1 by the late 2030s [63]. For the FCC-hh, defined by the target of 100 TeV proton-proton collisions, a total integrated luminosity of 20-30 ab −1 is considered [189].
In Tab. 16 the results are given for the final 4l signals. This time we consider higher H ±± mass of 1 TeV. The kinematic cuts are defined in section 5.4.
For pp collisions the 4e channel gives comparable to the background signals in MLRSM with v R = 15 TeV. In pp collision the lowest order t-channel is not present, so no destructive interference with the s-channel is possible. As given in Tab. 16, the maximal significance value is S = 11 [290] for 4µ in MLRSM with v R =15 TeV for HL-LHC and FCC-hh, respectively. For HTM in the same 4µ channel S < 1 both in NH and IH neutrino mass scenarios. So a detection of 4µ signals above the background level at HL-LHC and FCC-hh would give a clear indication for the MLRSM model with high values of v R .
So far, we have focused on comparisons of the two models looking for specific signals for the total four charged lepton production rates and compared it to the background processes. In this way, we can present clear differences in a prediction for the leptonic signals in both models. In particular, as shown in Tab. 15 and Tab. 16, there are cases where the SM background is comparable or exceeds the BSM signal for electrons and positrons in the finals state. The question is if in such cases, dilepton distributions can help to identify small BSM signals and discriminate further BSM models. In Fig. 14 the SM background and BSM dilepton distributions for e + e − → 4e are given. We consider distributions of pairs of electrons/positrons (the same charge leptons) and electronpositron pairs, assigned as SSDL and OSDL, respectively. As we can see, the SM invariant mass SSDL and OSDL distributions are quite uniform, in opposite to the BSM signals where clear peaks are present for SSDL signals. The reason is obviously that the same sign dileptons originate from the same doubly charged particle. As we can see, in this case, though 4e signals in Tab. 15 are below the SM background, both HTM and MLRSM dileptons can be identified. It is less visible for the leptonlepton separation ∆R ee though the SSDL (OSDL) signals are enhanced for higher (lower) values of ∆R ee , respectively, in both considered models. Let us note that m ee and ∆R ee distributions are very similar. This conclusion does not change for dimuon distributions or hadron colliders. For non-degenerate doubly charged masses in MLRSM in Fig. 14

Conclusions and outlook
The doubly charged Higgs bosons H ±± pair production at e + e − and pp colliders, with their subsequent decays to four charged leptons can give a very clear signal when searching for non-standard scalar particles effects without missing energy. We discuss a relation between vacuum expectation value of the triplet v ∆ and H ±± couplings with leptons, taking into account constraints on v ∆ coming from low energy studies connected with the ρ-parameter, muon (g − 2) µ , lepton flavor violation, e + e − , LHC processes, and neutrino oscillations (normal and inverse mass scenarios). The low energy experiments rule out v ∆ below 10 eV (for M H ±± ∼ 700 GeV) both for normal and inverted hierarchy, the strongest limit for non-zero mass of the lightest neutrino comes from LFV µ → 3e, see Still, assuming non-universality of leptonic decays, and due to fields richness of MLRSM, branching ratios for the H ±± decays can be very different in both models, leading to different final signals.
We discuss the same H ±± masses in both models. Taking into account all leptonic decays, we show that LHC experimental data still allow for H ±± mass as small as 700 GeV. We take it as the first scenario, the second is for H ±± mass equal to 1 TeV.
We discuss carefully possible decay channels and finally, we make predictions for the complete process pp → H ++ H −− → 4l. In both models, we optimised parameters to maximise separately e + e − (pp) → 4e and e + e − (pp) → 4µ signals, at the same time being in agreement with all experimental constraints coming from other considered processes.
The results are gathered in Tab. 15 and Tab. 16. There are many interesting conclusions that we can draw from them, as discussed in section 5.5. In general, due to kinematic cuts and chosen parameters, 4µ signals dominate over 4e. The latter signals are in most cases at best at the level of the SM background, both for lepton and hadron colliders. This situation gives a way to discriminate the two models. In fact, the most interesting situation in which v R in MLRSM is relatively large, above sensitivity of LHC (we took v R = 15 TeV) does not give too strict constraints on the model parameters, and the discovery signals can be large for e + e − (pp) → 4µ. In particular, for the HL-LHC and FCC-hh cases, detectable signals which would exceed the SM background are possible only for MLRSM. This conclusion is rather stable over changes of model parameters, for considered kinematic cuts. Though analysis of dilepton distributions can help further in detection of small BSM signals which are comparable or below the SM background, they are similar in patterns for both models and does not help in discrimination between HTM and MLRSM.
With the most straightforward setups, relying only on the production and decay total counting of events, we can discriminate models, and show in which channels we should look for that. We think that our work is an exemplary case study and from the minimal considerations, more sophisticated approaches can follow. As an outlook for further studies, a discussion of e ± e ± µ ∓ µ ∓ and e ± e ∓ µ ± µ ∓ channels might also be enjoyable, as well as four-lepton signal analysis with final state polarisation. It will be also interesting to investigate for chosen benchmark points processes with single produced H ±± or single charged Higgs scalars, and associated gauge bosons. For such cases the SM background will be much larger but it does not exclude positive BSM signals. 7 Appendix.

The HTM scalar potential and fields
The Higgs Triplet Model extends the Higgs sector of the SM by adding one scalar SU (2) L triplet (∆) with hypercharge Y = 2 to the Standard Model doublet Φ (following the convention Q = 1 2 Y + T 3 ). The most general scalar potential is given by [190] Without loss of generality we can take all the parameters to be real [191,192]. Denoting by v ∆ and v Φ the vacuum expectation values (VEV's) of the doublet and triplet We represent the scalar multiplets in the following way (7. 3) The triplet VEV v ∆ is expected to be at most at the order O (1) GeV to keep the electroweak ρ-parameter ∼ 1 [178,190,[193][194][195] (see section 3.2 for more details). The electroweak VEV is then given by The Yukawa sector contains the complete SM Yukawa Lagrangian along with an extra part for the triplet where, C is the charged conjugation operator, Y is the symmetric Yukawa matrix and are the left handed SU (2) doublets for the three lepton generations. After spontaneous symmetry breaking (SSB), the Yukawa couplings in Eq. (7.5) will lead to the Majorana mass matrix for the left handed neutrinos. The same term in the Lagrangian is responsible for the interaction between doubly charged scalar particles and charged leptons. The H ±± − l ∓ − l ∓ vertex breaks the lepton number (see sections 2,3). The fields, δ ±± = H ±± , represent the doubly charged scalar with the mass To get physical states for neutral and singly charged particles, appropriate rotation of fields in the CP-odd and CP-even sectors must follow Further, we use an approximation sin α ∼ 2 v∆ vΦ → 0 [86], neutral scalar masses becomes In the singly charged sector rotation of fields and masses are the following to obtain the charged Goldstone (G ± ) along with a singly charged scalar (H ± ) with mass 14) The H ± and H ±± scalar's squared masses (7.14) and (7.7) contain terms proportional to v 2 Φ and are inversely proportional to the triplet VEV v ∆ , which should be less than O (1 GeV) (see section 3.2). That means that M H ± , M H ±± can be at the level of a few hundred GeV or higher. Latest LHC bounds on the doubly charged scalar masses vary from 450 to 870 GeV, depending on the decay modes, assuming that BR (H ±± → l ± l ± ) ≥ 10% [55]. Photon-photon fusion studies [196] set a bound on M H ±± at the level of 748 GeV. Limits coming from e + e − colliders are significantly lower, from L3 Collaboration (LEP) it is about 100 GeV [102]. This bound comes with assumption that the t-channel is negligible (Fig. 1) as suppressed by the low H ±± − l − l coupling. For singly charged scalar masses the mass bound is even lower, M H ± = 80 GeV [124].
In this paper we assume that the neutral and charged scalars' masses are degenerated 3 , that means M H ±± = M H ± = M H = M A . That choice protects proper ranges of the T-parameter and potential unitarity for v ∆ 1 GeV [86,87].
where the quantum numbers in square brackets are given for SU (2) L , SU (2) R and U (1) B−L groups, respectively. The vacuum expectation values (VEVs) of the scalar fields can be recast in the following form: VEVs of the right-handed triplet (∆ R ) and the bi-doublet (φ), propel the respective symmetry breaking: The full scalar potential includes left and right-handed triplets [39,198,199]: ). (7.18) Though in HTM and MLRSM we have left-handed triplets, HTM is not a simple subset of MLRSM as the scalar potentials, SSB mechanism, VEVs and underlying physics which follows are different. The scalar potential (7.18) in MLRSM is much more complicated than its counterpart in HTM: in MLRSM the triplet ∆ L is intertwined with right-handed multiplet ∆ R and bidoublet φ. It makes relations among physical and unphysical Higgs boson fields rather complex in MLRSM. Here significant are relations between the α 3 scalar potential parameter (which includes a mixture of a bidoublet and triplet fields) and ρ 1 , ρ 3 scalar potential parameters for doubly charged Higgs boson masses given in Eq. (7.28) and Eq. (7.29) below. In correlation with experimental constraints for singly charged and neutral scalar fields, these parameters give the lowest limits for doubly charged Higgs masses, as discussed in [169]. Moreover, due to Yukawa couplings of left-and right-handed leptons with bidoublet in Eq. (7.38), the doubly charged Higgs bosons in both models couple differently to leptons. Consequently, both model neutrino mass relations are different, in HTM restricted directly by neutrino oscillation data, as discussed in the main text.
After spontaneouss symmetry breaking of the potential Eq. (7.18), the mass matrix which includes M H 0 0 can be written in the following form (for details, see [198]) Expanding eigenvalues of this matrix in a small = κ 2 The analytic mass formulas for other scalar bosons in MLRSM as a function of quartic couplings and v R can be written as [56] M 2 where κ 1 , κ 2 are VEVs of the bidoublet and κ 2 1 + κ 2 2 has to be equal to the electroweak symmetry breaking scale v, see (7.4). We assume that κ 1 = v = 246 GeV and κ 2 → 0. Some explicit masses of Higgs bosons relevant for H ±± branching ratios in section 5.3 comes from restrictions discussed in [56].
In MLRSM relations among physical and unphysical fields ("G" stands for Goldstone modes) are The structure of Higgs potential in a general framework of left-right symmetric models has been discussed in details in [198,199]. We adopted it in our studies. In particular, to retain the invariant Majorana Yukawa couplings of the leptons to the Higgs triplet, the potential does not include the terms with all multiplets (bidoublet, two triplets) present simultaneously, e.g. [198,199] denoted as the β i -type terms). In the limit of vanishing β i terms the doubly charged Higgs scalar 2 × 2 mass matrix is diagonal and does not include the mixed mass terms δ ±± L δ ∓∓ R . These restrictions simplify a form of doubly charged mass terms, as given in the manuscript, Eqs. (7.28) and (7.29). It means that in MLRSM (and other extensions when gauge couplings g L = g R ), doubly charged Higgs triplets in ∆ L,R are physical fields. There is no mixing angle between two doubly charged Higgs bosons in MLRSM and the mass matrix which appears there for unphysical fields is diagonal from the very beginning. This no-mixing feature is also true in general where the β i -type terms are allowed, in the limit v R κ 1,2 . In MLRSM due to additional heavy states the neutrino sector and Yukawa couplings are more complicated than in HTM. Here we argue that due an energy scales difference between v R and the low-energy bidoublet VEV κ ≡ κ 2 1 + κ 2 2 , κ v R , the see-saw mechanism is possible and low energy LFV signals are suppressed due to high v R and heavy neutrino masses. To see, this, the most general doubly charged couplings to leptons, which takes into account mixing matrices, reads [39] δ ++ . These couplings originate from the Yukawa part of the Lagrangian for additional scalar triplets and a bidoublet φ: Uniqueness of left-and right-handed couplings for positively and negatively charged doubly charged Higgs bosons to leptons in Eq. (7.36) and Eq. (7.37) is due to the Feynman rules and flow of the charged currents in vertices, as explained in length in [200]. The relations between physical and unphysical scalar, gauge and fermion fields are embedded in our FeynRules package to calculate branching ratios and cross sections. We also mentioned that for neutral scalars, due to the bidoublet coupling in Eq. (7.38) to both left-and right-handed leptons, there is a mixture of scalar fields coming from left and right triplets, however, as given in Eqs. (7.30)-(7.35), for v R κ 1,2 , most of the mixings are negligible.
Diagonalization of the resulting neutrino mass matrix goes with the help of a unitary 6×6 matrix U This procedure leads to the introduction of the K L and K R submatrices in Eq. (7.36) and Eq. (7.37) [145,201]. The charged lepton mass matrix is diagonalised by V l Apart from charged lepton and neutrino mass terms, Lagrangian Eq. (7.38) contains scalar-lepton interactions too. (M ν ) diag contains 3 light neutrinos, there contribution to the couplings Eq. (7.37) and Eq. (7.36) are negligible. To see amount of heavy neutrinos contributions to Eq. (7.37) and Eq. (7.36), we note that structure of K L and K R mixing matrices are the following [145,201] (K L ) liνj =            e µ τ   · · · · · · · · ·      light neutrinos   · · · · · · · · ·      heavy neutrinos Off-diagonal elements for heavy neutrinos couplings in K R are typically also of the order of inverse heavy neutrino mass scale, that is why LFV couplings of leptons with doubly charged Higgs bosons are strongly suppressed, δ L,R − l − l off-diagonal lepton couplings are suppressed by 1/m 2 N when comparing to diagonal cases.
For reasons discussed in [149] and more extensively in [202], we take seesaw diagonal light-heavy neutrino mixings. It means that W 1 couples mainly to light neutrinos, while W 2 couples to the heavy ones.
To summarize, unlike in the HTM case, the H ±± − l ∓ − l ∓ vertex does not depend on the light neutrino mixing. With v L = 0 the MLRSM realizes the seesaw type-I mechanism and the light neutrino mass is due to the existence of additional heavy neutrino states and v R scale.
We should note that it is not natural and very hard to create non-decoupling mixings for nondiagonal K L and K R matrix elements, even when some symmetries are considered in type-I seesaw models [202].

Supplemental material for phenomenological studies of H ±± scalar particles
Diagrams in Fig. 4 present the contribution from the singly and doubly charged scalar particles to lepton flavour violating processes and to muon (g−2) µ . Those diagrams contain vertices H ±± −l i −l j and H ± − l i − ν j which origin from the Yukawa part of the Lagrangian, combining the Standard Model Yukawa term with the triplet part Eq. (7.5) From Eq. (7.42) we obtain the interaction between charged leptons and a doubly charged scalar and the interaction of a singly charged scalar with a charged lepton and neutrino. Taking into account Eq. (7.13) and Eq. (7.4) and keeping in mind that y ij ∝ 1 vΦ is a SM diagonal matrix Vertices V ± ∆ and V ±± ∆ comes from the same part of the Lagrangian and they break the lepton flavour. Vertex V ± Φ is proportional to v ∆ while vertex V ± ∆ is inversely proportional to the triplet VEV and dominates up to v ∆ ∼ 10 6 eV. Since we are interested in lower regions of v ∆ values, its effect is negligible. So, with a good approximation for low values of v ∆ As we discussed in section 3, the branching ratios of the radiative and µ-to-e conversion depends on the one-loop form factors. From Eq. (3.5), we can read them explicitly We have checked with earlier literature and the analytic forms for the CLFV processes are given as [33,92,93,105,203]: Radiative lepton decay l i → l j γ: The branching ratios of radiative decay processes can be given by: BR(l i → l j γ) = 384π 2 (4πα em )|A R | 2 BR(l i → l j ν li ν lj ).
The contribution of H ±± to the branching ratios is eight times larger than by H ± because of the difference in a magnitude of couplings between V ± and V ±± in Eq. (7.46), in addition the amplitude is proportional to the particles charge (which gives an additional factor of 4).
Three body decays l → l i l j l k : In the computation of conversion rate µ to e, both form factors contribute and the analytic form can be written as (7.50) Therefore, the µ-to-e conversion ratio in the nuclei field can be given as [93,105]: Al 0.7054 × 10 6 0.4643 × 10 −9 11.5 Table 17: Total muon capture rate and effective charge for 197 Au, 48 Ti and 27 Al [204]. .

Muon (g − 2) µ
In case of doubly charged scalars, cumulative effects of muon (g − 2) µ and lepton flavor violation have been discussed in details in [105,205,206] and for a triplet scalar it has been discussed in [91].
f l is a symmetric factor equals to 4 for l = µ and 1 otherwise. The term proportional to the C 1 integral is connected with the diagram a) in Fig. 15, the term with the C 2 integral corresponds to Fig. 15 b). Equation (7.54) presents a contribution from a singly charged particle H ± (Fig. 15 c) ) to (g − 2) µ . Since both V ± and V ±± vertices are comparable, see Eq. (7.46), the contributions from different diagrams depend mostly on C 1 , C 2 and C 3 integrals. Fig. 16 shows that the strongest contribution (g − 2) µ comes from the doubly charged scalar H ±± .