Probing Zee-Babu states at Muon Colliders

The Zee-Babu model is a minimal realization of radiative neutrino mass generation mechanism at the two-loop level. We study the phenomenology of this model at future multi-TeV muon colliders. After imposing all theoretical and low-energy experimental constraints on the model parameters, we find that the Zee-Babu states are expected not to reside below the TeV scale, making it challenging to probe them at the LHC. We first analyze the production rates for various channels, including multi singly-charged and/or doubly-charged scalars at muon colliders. For concreteness, we study several benchmark points that satisfy neutrino oscillation data and other constraints and find that most channels have large production rates. We then analyze the discovery reach of the model using two specific channels: the pair production of singly- and doubly-charged scalars. For the phenomenologically viable scenarios considered in this study, charged scalars with masses up to ${\cal O}(3$--$4)$ TeV can be probed for the center-of-mass energy of $10$ TeV and total luminosity of $10~{\rm ab}^{-1}$.


I. INTRODUCTION
Although the Standard Model (SM) of particle physics is the most successful theory to describe nature at the fundamental scale, it has several drawbacks.Among its shortcoming, the most prominent is that neutrinos are massless in the SM.However, several experiments [1][2][3][4][5][6][7] have discovered neutrino oscillations that firmly established non-zero neutrino masses for at least two generations.Observations of non-zero neutrino masses unquestionably call for physics beyond the SM (BSM).A well-motivated class of models explaining the origin and the smallness of the neutrino mass is the radiative neutrino mass [8] generation mechanism where the neutrino mass is generated at one or more loops with the BSM particles whose masses are typically not too above the TeV scale.An economical way of generating neutrino masses, which only extends the scalar sector of the SM, are the Zee model [9,10] and the Zee-Babu model [11,12].In the former case, a singly-charged scalar and a second Higgs doublet are introduced, and neutrino masses appear at the one-loop level.The latter model extends the SM scalar sector by only a singly-charged scalar and a doubly-charged scalar, and, consequently, neutrino masses arise at two-loop order.This model has been widely studied in the literature [13][14][15][16][17][18][19] (see also Ref. [20] for a comprehensive review), and different versions are also proposed.For example, colored versions of the Zee-Babu model are studied in Refs.[21][22][23][24][25][26].
In this work, we focus on the minimal version of the Zee-Babu model and consider the possibilities of observing the new physics states, namely, the single-charged and doubly-charged scalars at the muon colliders.A salient feature of this model is that owing to two-loop suppression, obtaining the correct neutrino mass scale requires that the new physics states are not too heavy or the Yukawa couplings are not too small.For the BSM scalars of masses O(100) GeV, there is a lower bound of > 10 −2 on the largest of the Yukawa couplings.Consequently, several lepton flavor violating (LFV) processes mediated by the Zee-Babu scalars can become observable.Among them, the most promising LFV process is µ → eγ, which cannot be arbitrarily small.As will be explained later, if fine-tuned regions of the parameter space are not considered, then non-observation of LFV signals and reproducing the observed neutrino oscillation data prefer the new physics states to reside at or above ∼ TeV.Moreover, a recent analysis of the LHC puts a lower limit of about TeV on the mass of the doubly-charged scalars.
Recently, multi-TeV muon colliders have attracted interest in discovering new physics at the TeV scale and beyond (for a review see Refs.[27][28][29][30]).One of the important aspects of multi-TeV muon colliders is that we can achieve both a high center-of-mass energy and a clean environment since the backgrounds are mostly of electroweak origin.For example, a signal-to-background ratio for Higgs boson production at muon colliders is about 10 −2 for all the center-of-mass energies while it is about 10 −6 at the LHC.This would imply that the multi-TeV muon colliders are the perfect colliders for discovery, precision measurements, and even new physics characterisation.Phenomenological analyses have been extensively carried out at future muon colliders both for the SM and new physics .In this work, we show that multi-TeV muon colliders have great potential to discover charged scalar states.One of the main reasons for this is that new physics states typically couple with muons stronger than the other generations within this framework.Moreover, the doubly-charged scalar may leave some clean signatures through which this model can be efficiently probed at future muon colliders.In particular, using a fully-fledged analysis at the detector level, we show that singlyand doubly-charged scalars of masses about 1. 25-3.75 TeV can be efficiently probed at muon colliders.
The study performed in this work for searching Zee-Babu states in muon colliders may also apply to various models.Because, apart from the Zee and Zee-Babu models, there is a plethora of radiative neutrino mass models for which charged scalars are essential ingredients.The scotogenic model [56] that links dark matter to the radiative generation of neutrino masses comes with a singly-charged scalar originating from an inert doublet.Moreover, going beyond two-loop, there are three minimal three-loop neutrino mass models, namely (i) KNT model [57], (ii) AKS model [58], and (iii) cocktail model [59].The KNT model consists of two singly-charged scalars and three copies of right-handed neutrinos.The AKS model extends the SM scalar sector with a second Higgs doublet, a SM singlet, and a singly-charged scalar.In addition to the Zee-Babu states, the cocktail model introduces an inert doublet.In this model, however, the singly-charged scalar carries a discrete charge, and the Zee-Babu loop is not allowed.Another model worth mentioning is the BNT model [60], which, in addition to a singly-charged and a doubly-charged scalars also predict the existence of a triply-charged scalar.All the models mentioned so far assume neutrinos are Majorana particles; however, neutrinos can also be Dirac.The minimal radiative Dirac neutrino mass models also predict the existence of charged scalars [61][62][63][64]).The above list does not provide a complete set of references, and we refer the reader to the recent review Ref. [65].
This work is organized in the following way.After introducing the model in Sec.II, we discuss various theoretical and experimental constraints on model parameters in Sec.III.By imposing the derived constraints, we perform numerical fit to the neutrino oscillation data in Sec.IV and present five benchmark case studies.A detailed analysis of the production rates of new scalar states at muon colliders is then carried out in Sec.V and discovery prospects for two production channels are illustrated in Sec.VI.Finally, we conclude in Sec.VII.

II. MODEL
The details of Zee-Babu model [11,12] are presented in this section.The particle content of the SM is extended with only a singly-charged and a doubly-charged scalars, The kinetic part of the Lagrangian for these two fields is given by, The scalar potential of this theory also takes a simple form, where H ∼ (1, 2, 1/2) denotes the SM Higgs.The last term in V NP plays a crucial role in generating Majorana neutrino masses.Moreover, the Yukawa couplings associated to these new states take the following form: Here, we have used the notation L ∼ (1, 2, −1/2) ∼ (ν L , ℓ L ) T and ℓ ∼ (1, 1, −1) = ℓ R , and i, j (a, b) are family (weak) indices.The 3×3 Yukawa coupling matrices f and g are anti-symmetric and symmetric, respectively, in the family space.Therefore, f has three independent entries (f eµ , f eτ , f µτ ) and g has six (g ee , g eµ , g eτ , g µµ , g µτ , g τ τ ).For the simplicity of the analysis, we take all parameters to be real.
Neutrino mass generation at the two-loop order.
If we assign two units of lepton number (L) to the BSM states, i.e., L[ϕ] = −2, L[κ] = −2, then all terms in the entire Lagrangian conserve lepton number except the term ϕ + ϕ + κ −− in the scalar potential.Consequently, lepton number is broken by two units and neutrinos acquire Majorana mass at the two-loop level as depicted in figure 1.
This two-loop diagram is calculable, and neutrino mass formula takes the following form [15]: where charged lepton masses are denoted by m i .Since the masses of the charged leptons are much smaller than the charged scalar masses, they can be neglected from the corresponding propagators in the loop integral.Then the corresponding loop function becomes essentially generation independent, and is given by, with In some special cases, the loop function takes the following simple forms: For numerical analysis, we evaluate the integral given in equation ( 8), and use the following expression for the neutrino mass matrix: In the above equation, m diag E represents the diagonal charged lepton mass matrix, namely m diag E = diag(m e , m µ , m τ ).Note that due to the antisymmetric nature of the Yukawa matrix f , the determinant of the neutrino mass matrix vanishes.Therefore, lightest of the neutrinos remain massless at twoloop order.In general, the neutrino mass matrix given in equation ( 6) can satisfy both the normal and inverted mass ordering.Neutrino oscillation data for these two cases are collected in Table I.Interestingly, due to the two-loop suppression, the larger of the couplings have to be ≳ 10 −2 to reproduce the correct neutrino mass scale.Consequently, this model can be efficiently probed through the rare charged lepton flavor violating processes.For analysing the implications in muon colliders, we are particularly interested in coupling of order ∼ O(1), especially the µµ coupling to new physics states.Before we perform a fit to the neutrino masses and mixings, in the next section, we summarize the theoretical and experimental constraints on the model parameters.

III. THEORETICAL AND EXPERIMENTAL CONSTRAINTS
In this section, we summarize various theoretical as well as experimental constraints on the model parameters.

A. Theoretical constraints
Model parameters that play role in neutrino mass generation are the Yukawa couplings f ij , g ij and the cubic coupling µ.Therefore, we are interested in finding the constraints on these parameters.First, due to the perturbativity of the Yukawa couplings, we restrict ourselves to f ij , g ij ≤ 1.As for the cubic coupling, the most crucial bound arises from loop effects.In  particular, due to the cubic coupling ϕ + ϕ + κ −− , effective quartic interactions for ϕ + and κ ++ are generated at one-loop level.These effective quartic coupling are negative, and with m ϕ , m κ masses of similar order, they are given by λ eff ∼ −µ 4 /(6π 2 m 4 ϕ ) [13].Hence, the tree-level quartic couplings (λ ϕ , λ κ , λ ϕκ ) must be positive and larger in magnitude such that the net effective quartic couplings are positive to ensure the vacuum stability.Then, assuming perturbitivity and restricting to values λ ϕ , λ κ , λ ϕκ ≤ 1, constraints on the cubic coupling parameter is obtained [13], which read Since all the cases we study are such that m κ and m ϕ are closer in mass, while performing a numerical analysis, we impose the first constraint listed above.
As will be shown in the next section, a fit to the neutrino oscillation data prefers somewhat larger values of µ.Such a large value compared to the masses of the scalars can be dangerous if it leads to a deeper minimum of the scalar potential for non-vanishing values of charged scalars leading to charge-breaking minimum.One must impose relevant constraints to avoid charge breaking minimum, which, for example, are studied in Refs.[67][68][69].By applying similar method to Zee-Babu model, a conservative bound on the cubic coupling is obtained in Ref. [18], However, a much more restrictive bound on this parameter is suggested in Ref. [18] (see also Ref. [17]), which is given by, In our numerical analysis, we make sure that all these conditions, equation (11) (the first line), equation (12), and equation ( 13) are met.[73,74].

B. Low-energy experimental constraints
The Yukawa couplings that give rise to neutrino masses also participate in rare charged lepton violating (cLFV) processes.Such cLFV processes ℓ − i → ℓ + j ℓ − k ℓ − l already take place at tree-level and ℓ − i → ℓ + j γ appears at one-loop level.The partial widths for these decays are given by [15], and, In the last formula, δ kl takes into account the possibility of having two identical particles in the final state.Utilizing the above formulas, the relevant branching fractions are obtained by computing BR(ℓ At one-loop, new physics contribution to the anomalous magnetic moment of the i-th generation of lepton reads [15], Furthermore, another LFV process, namely, the muonium-antimuonium transition is also mediated via the tree-level interactions.Consequently, the following effective Lagrangian coupling coefficient is generated [15,72] The interactions of the singly-charged scalar generate four-fermion operators that lead to a charged current involving four leptons ℓ i → ℓ j νν.These processes are highly constrained from lepton universality tests.The current constraints are given by [75], We impose 2σ constraints on these quantities, and δ is defined as, where as usual v = 246 GeV.Furthermore, not to change the Fermi constant from the SM prediction, we impose 2δ(µ → eνν) < 0.002 [76].
The singly-and the doubly-charged scalars do not interact with the quarks; therefore, lepton flavor violating µ − e conversion in the nuclei is forbidden at the tree-level.Such processes are then induced at one-loop level and are strongly correlated with µ → eγ.However, as long as the latter mode satisfies the current experimental bounds, the µ − e conversion is also below the present bound.This is why we do not include the corresponding bounds in our analysis.
The current experimental bounds and the associated constraints on the model parameters arising from tree-level mediated cLFV processes are summarized in Table II.Similarly, the current limits on cLFV decays and lepton anomalous magnetic moment that happen at the one-loop level along with the constraints are recapitulate in Table III.

C. Current collider constraints
At the LHC, both the singly-and the doubly-charged scalars are pair produced via Drell-Yan mechanism.Recently, ATLAS collaboration has presented their updated search on the doubly-charged scalar in the same-charge two-lepton invariant mass spectrum.In the context of Zee-Babu model, their study includes e ± e ± , e ± µ ± , and µ ± µ ± final states with the identification of three/four leptons.With 139 fb −1 of data from proton-proton collisions at √ s= 13 TeV, no significant excess above the SM prediction was found.This rules out the Zee-Babu state, κ ±± , with masses smaller than 900 GeV [77], whereas, the projected sensitivity with 3 ab −1 of data at √ s= 13 TeV, the limit is 1110 GeV [78].However, there are no dedicated searches for the weak-singlet singly-charged scalar at the LHC.Recently, Ref. [79] performed an analysis of such scalars and found that a charged scalar decaying exclusively into e ± and µ ± can be excluded up to masses of 500 GeV with an integrated luminosity of 200 fb −1 data.

IV. NUMERICAL ANALYSIS
In this section, we perform a numerical fit to the neutrino oscillation data.The neutrino mass matrix is given in equation ( 6) that contains nine Yukawa couplings f ij , g ij and a cubic coupling constant µ.Without loss of generality, the charged lepton mass matrix can be made real, diagonal, and positive.In this basis, via field redefinition, one can make elements of f real.By further field redefinition, one of the phases from g can be removed as well as the cubic coupling can be made real.This leaves us with ten real parameters and five phases.However, for the simplicity of our analysis, we consider all parameters to be real (i.e., set all phases to zero).For the masses of the charged leptons appearing in the neutrino mass formula, we use the PDG values [80].In the fitting procedure, we make sure that the parameters (f ij , g ij , µ) satisfy all the theoretical constraints mentioned above.Furthermore, the loop function entering in neutrino mass formula involves masses of the two scalars, which we fix to certain values.
This numerical study is based on a χ 2 analysis, and which is defined as, Here σ i represents experimental uncertainty; E i and T i stand for experimental central value and theoretical prediction for the i-th observable respectively.In the above equation, i is summed over five observables, two neutrino mass-squared differences and three mixing angles.Experimental values of these quantities are given in Table I.In this fitting procedure, we do not include the CP violating phase, δ, in the χ 2 -function.This phase has not been measured yet in the experiments, and a global fit to the neutrino oscillation data currently allows almost the entire range from 0 to 2π.Since the parameters, for a fixed set of masses (m ϕ , m κ ), are severely constrained by the flavor violating processes, we perform a constrained minimization.In this procedure, the χ 2 function is subject to all constraints, theoretical as well as experimental, discussed in Sec.III.Since it is highly challenging to explore the entire parameter space, we consider five benchmark points.Fit parameters obtained for the benchmark points are listed in Table IV.Moreover, in the same table, we summarize various branching ratios for ϕ ± and κ ±± decay modes.The outcome of these benchmark fits are presented in Table V.In this work, we only focus on NO for the neutrino masses, which demands the largest Yukawa coupling to be g µµ , hence has profound implication in the muon colliders.If instead, IO is considered, the largest Yukawa coupling becomes g µτ (see, for example, Ref. [18]).
From our numerical analysis, we find that the most stringent constraint on the Yukawa parameters and the masses of the Zee-Babu states arise from the charged lepton flavor violating µ → eγ process.In particular, our result shows that without going to a fine-tuned part of the parameter space, the minimum masses possible for the Zee-Babu states is close to ∼900 GeV, which is fully consistent with the previous findings [17,18].Fine-tuned regions of the parameter space do allow Zee-Babu states to have masses as low as ∼ O(200) GeV [18].It is important to note that even if other important flavor violating processes, for example, µ → 3e, may be suppressed by taking vanishing g ee or g eµ , however, µ → eγ cannot be suppressed to arbitrarily small values since both the Yukawa couplings f and g contribute to it that cannot taken to be zero simultaneously.We emphasise that it is still difficult to fully rule out the Zee-Babu model from cLFV constraints.Moreover, the BSM scalars can be as heavy as ∼ O(100) TeV   Here we show also some of its characteristics like the decay branching ratios of the singly-charged (ϕ ± ) and doubly-charged (κ ±± ) scalars.These benchmark parameters along with the neutrino mass matrix can be provided as a text file upon request.
depending on perturbativity conditions imposed; for a detailed discussion on these, we refer the readers to Refs.[17][18][19].
A fit to the neutrino oscillation data requires the largest entry to be the µµ coupling with the doubly-charged scalar.The coupling of the singlet is also expected to be sizable with the muon.This is why, muon colliders are the perfect machines to search for these states that  have masses in the multi TeV range.In our studied five benchmark scenarios, we have fixed the masses of these states to be (m ϕ , m κ ) = (m, m), (2m, m), (m, 2m), (3m, m), and (m, 3m) (with m = 1250 GeV), which have great potential to be discovered in the muon colliders.All these benchmark scenarios also predict large flavor violating processes in the muon and tau decays, specifically µ → eγ, µ → 3e and τ → 3µ, τ → e + µµ modes are lying just below the current limit and will be observed in the near future, see Table V.The future sensitivities of these flavor violating processes are BR∼ 6 × 10 −14 for µ → eγ [81], BR∼ 10 −16 for µ → 3e [82], and BR∼ 10 −9 for τ → 3µ, τ → e + µµ [83].These promising flavor violating processes that are within reach of near future experiments are highlighted with boldface in Table V.
Before concluding this section, we remark that another interesting neutrino mass model, namely, the type II seesaw that leads to neutrino mass at tree-level, also predicts a singlyand a doubly-charged scalars arising from a scalar weak-triplet.Phenomenology of this model is different from that of Zee-Babu model and there are much more freedom in choosing the Yukawa couplings to fit oscillation data.A comprehensive study of type II seesaw model and distinguishing between these two models at muon colliders will be presented in a future work.

V. THE ZEE-BABU MODEL AT MUON COLLIDERS: PRODUCTION RATES
In this section, we analyse the production rates for the Zee-Babu states at multi-TeV muon colliders.The different production channels can be categorised into three classes: (i) the production of charged lepton pairs, (ii) the production of singly-charged scalars, and (iii) the production of doubly-charged scalars.For all the calculations of the production rates, we make use of Madgraph5 aMC@NLO version 3.4.2[84] with the use of the publicly available FeynRules [85] model file which can be found in this link https://feynrules.irmp.ucl.ac.be/wiki/ZeeBabu.

A. Charged lepton pair production
For charged lepton pair production, we consider both the flavour-conserving (FC) as well as the flavour-violating (FV) channels (see figure 2).The FC lepton pair production cross section is shown in figure 3. The contribution of the doubly-charged scalar to the production cross section for the lepton flavour ℓ is proportional to g 2 ℓℓ .As per the benchmark points we have chosen a noticeable effect can be seen in the muon channel since the g κµµ coupling can be of order O(1).While a suppression of the muonic cross section is caused by the ∝ m ℓ /m 2 κ factor.Such effects can be seen clearly in the right panel of figure 3. We can see that the effects of the doubly-charged scalars are mainly driven by the value of the g µµ coupling since we can see that the predictions for m κ = 2.5 TeV (green) and m κ = 3.75 TeV (blue) are exactly the same because they correspond to the benchmark points BP3 and BP5 for which g µµ ≈ 1.For the case of the electron channel, the contribution of the doubly-charged scalar is negligibly small since the corresponding coupling is of order O(10 −2 -10 −1 ) and therefore the value of the corresponding cross section is mainly driven by the contribution of the gauge-boson Feynman diagrams.The contribution to the τ -lepton pair production is even much more smaller for which reason we do not show the corresponding result in this paper.The cross section for the FV lepton pair production is shown in figure 4. Contrarily to the FC case, the cross section in this case receives only contribution from the t-channel diagrams and has the following behavior: From the choice of the couplings g ij in our benchmark scenarios, it is expected to have σ µ ± τ ∓ ≫ σ e ± µ ∓ ≫ σ e ± τ ∓ a feature that can be clearly seen in figure 4. We note that the cross section for the e ± τ ∓ channel is orders of magnitudes smaller and therefore we dot not show it here.We note that these processes are strongly correlated to charged lepton flavour violating decays.However, for our benchmark scenarios the most important probe would be the search for the production of µ ∓ τ ∓ given that the other two are suppressed and therefore we do not expect them to give additional information.(red), BP2 (green), BP3 (blue), BP4 (cyan) and BP5 (orange).Note that for the production of the e ± µ ∓ , the cross section in BP1 is not shown as it is well below 10 −8 fb.

B. Production of singly-charged scalars
We turn now into the discussion of the production mechanisms of the singly-charged scalars at muon colliders.We analyse two production channels: (i) the pair production mode (upper panel of figure 5) and (ii) the production of four scalars (lower panel of figure 5).The production rate for singly-charged scalar pairs receive two contributions: a Drell-Yan type contribution (left upper panel of figure 5) and t-channel contribution (right upper panel of figure 5).We note that the second contribution is always suppressed as was checked out since the relevant couplings -f eµ and f τ µ -are small in all the benchmark scenarios under consideration.The production cross section for ϕ + ϕ − fall as 1/s ( √ s is the center-of-mass energy) as can clearly be seen in the left panel of figure 6 with the maximum being about 0.5 fb for the benchmark points corresponding to m ϕ = 1.25 TeV (BP1, BP3 and BP5).We can see also that the total cross section for this mass value does not depend on the values of the Yukawa-type couplings which is another test of the negligible contribution of the t-channel diagrams in this process.
For the other benchmark points, the production cross section reaches a maximum just above the threshold and then decreases as 1/s with the cross section for BP2 (m ϕ = 2.5 TeV) is always larger than for BP4 (m ϕ = 3.75 TeV).
For the production of four singly-charged scalars, there are multiple contributions that can be categorised into four classes: (i) the gauge-boson contribution through the exchange of a photon or a Z-boson and which occurs always in s-channel, (ii) the contribution of charged leptons and neutrinos through t-channel which is shown on the left lower panel of figure 5, (iii) the contribution of the SM Higgs boson with a neutrino exchanged in the t-channel and (iv) the contribution of a doubly-charged scalar.We note that the contribution of the leptons is always subleading since it involves at three propagators and a factor of f 2 µℓα f 2 µℓ β f 2 ℓαℓ β which implies a suppression of about O(10 −8 -10 −12 ).The second contribution of the SM Higgs is proportional to the scalar quartic coupling λ ϕ which we have chosen to be equal to one.The last contribution is the most dominant as it is proportional to µ 4 and which has been chosen to be large in this analysis.Note that in this case, the cross section get threshold enhancement in scenarios where m κ ≥ 2 m ϕ which is the case for BP3 and BP5.We can see clearly these features in the right panel of figure 6 as the cross section for these benchmark points can reach up to 30 fb near the production threshold.For BP1 where m κ = m ϕ the cross section is about three orders of magnitude smaller.For heavy singly-charged scalars (BP2 and BP4) the cross section do not go above ≈ 10 −4 fb.

C. Production of doubly-charged scalars
We now analyse the production of doubly-charged scalars at muon colliders.Similarly to the production of singly-charged scalars we study both the production of two scalars and of four scalars.Examples of the Feynman diagrams are depicted in figure 7. The pair production of doubly-charged scalars proceeds through the contribution of Drell-Yan s-channel diagrams (left upper panel of figure 7) and through t-channel diagram with the exchange of charged leptons (right upper panel of figure 7).Unlike the singly-charged case, the contribution of the t-channel diagram to the total rate of doubly-charged pair production is very important thanks to the magnitude of the g µµ coupling.In fact, we have checked that this contribution is the FIG. 8. Same as in figure 4 but for the production of κ ++ κ −− (left) and of κ ++ κ −− κ ++ κ −− (right).
most important for all of the benchmark points; for example, it is about 20 times larger than the s-channel Drell-Yan for BP1 at √ s = 3 TeV and becomes much larger for higher energies.
On the other hand, the t-channel contribution is found to interfere destructively with the schannel contribution (I ≡ 2Re(M * t M DY )/|M total | 2 ≈ −0.24 at √ s = 3 TeV for BP1).In the left panel of figure 8 we show the rate of the pair production cross section of doubly-charged scalars as function of √ s.We can see that the production cross section in this case is more important for BP3 and BP5 for which case g µµ is equal to 1.It falls with the center-of-mass energy but with a minimum of about ≈ 2 fb for BP4 which is a clear sign of the effect of the t-channel diagrams.
We close this section by a brief discussion of the production rate of four doubly-charged scalars at muon colliders.The Feynman diagrams are displayed on the lower panel of figure 7 and the production cross section is shown on the right panel of figure 8.In addition to the gauge-boson contribution, this process receives similar contribution as in the case of the singly-charged scalars except that there is no contribution from the µ-term.This would imply that the total cross section is smaller than in the case of four singly-charged scalars.We can see from figure 8 that the corresponding cross section varies in the range of 10 −3 -10 −1 fb.Two important remarks are in order here.First it seems from the same figure that the cross section violates unitarity but we have checked that basically for energies above √ s ≈ 200 TeV the cross section starts to decrease.Finally, despite the smallness of the corresponding cross section this channel leads to a final state of eight highly-energetic muons, which can be considered as golden channel for discovering the doubly-charged scalar at muon colliders.

VI. SENSITIVITY REACH
In this section we discuss the prospects of detecting Zee-Babu states at multi-TeV muon colliders.We will perform a simple estimates of the signal-to-background ratios at both the parton and the detector levels.For this task, we will only study two channels: the pair production of the singly-charged scalars and the doubly-charged scalars.In this study we consider the sensitivity reach for √ s = 10 TeV and L = 10 fb −1 assuming that the integrated luminosity increases linearly with the center-of-mass energy to compensate for the 1/s decrease in the cross section of the most processes that we can encounter at future muon colliders.

A. Technical setup
We consider the sensitivity reach of the Zee-Babu model at future muon colliders in the ϕ + ϕ − and κ ++ κ −− channels.For each case, we consider first a simple analysis at the parton level by comparing the signal and the background cross sections as a function of the minimum cut on the lepton transverse momentum.To get realistic estimates of the expected sensitivity reach, we also perform a fully-fledged signal-to-background optimization taking into account all the parton-shower and detector effects.In the detector-level part of the study, we consider only the muonic decay channels, i.e.
In this case, we have signatures consisting of two muons and missing energy for the ϕ + ϕ − channel and four muons for the κ ++ κ −− channel.In this subsection, we describe in detail the technical setup, including the Monte Carlo event generation and the detector simulation setup.For the background processes, the main contribution arises from the production of two massive gauge bosons.The 2 muons plus missing energy signature receives background contributions from the production of W + W − , W ± Z, and ZZ.For all these backgrounds, we consider both the production through the µ + µ − annihilation and the VBF fusion, which is estimated using the prescription of Ref. [51].The cross sections for the W W/W Z and ZZ processes, including the decays into muons and neutrinos, are 0.65 fb and 0.06 fb, respectively.In the case of four muon signature relevant to the pair production of doubly-charged scalars, the only background contribution arises from the production of ZZ, which can be categorized into two categories: s-channel like production of ZZ and VBF production of ZZ with cross sections of order 10 −3 -10 −1 fb including their decay branching ratios into muons 1 .Events are generated using Madgraph aMC@NLO version 3.4.2[84] and then passed to Pythia version 8310 [86] to add parton showering and hadronization.To reduce the effects of statistical fluctuations in the tails of the kinematical distributions, we have generated about 3-5 million events for the background processes depending on the respective cross sections.To take into account detector effects, we use the Simplified Fast Simulator (SFS) module [87] in MadAnalysis 5 [88][89][90][91][92][93].
The smearing and identification efficiencies of charged leptons and photons are implemented using the detector projection shipped with Delphes version 3.4.0[94] 2 .For charged muons, we require them to be tightly isolated in the sense that the sum of the transverse momenta of all the tracks within ∆R = 0.2 of the muon candidate excluding the muon itself to satisfy: To estimate the sensitivity reach, we compute the signal significance defined using the Asimov formula [96] On the other hand, we also consider the case where the background yields have some uncertainty denoted by δ.In that case, the significance formula in equation 26 becomes where n s and n b refer to the number of signal and background events, respectively, and δ b = xn b is the uncertainty on the background yield.
B. Sensitivity reach in the ϕ + ϕ − channel

Parton level
We start by the pair production of the singly-charged scalars at future muon colliders.From the choices of the parameters in the different benchmark points, we found that the outcome of BP1, BP3 and BP5 is the same since they correspond to m ϕ = 1.25 TeV.The corresponding branching ratios of the charged scalars are quite similar in these three benchmark points (see table IV).Therefore, we will study the sensitivity reach for BP1, BP2, and BP4 at √ s = 10 TeV.At the parton level, we consider three final states: The dominant backgrounds to these production channels are the diboson production: W + W − and ZZ where in the first background, both the two gauge bosons decay leptonically while in the second, only one of them decays leptonically and the other decays invisibly.The number of signal and background events is defined for the µ ± e ∓ + E miss T as    where the factor 2 is included to take into account combinatorics.Similarly, we can define the number of signal and background events for the case of the µ ± τ ∓ + E miss T channel as and of the µ ± µ ∓ + E miss T channel as where in the calculation of the number of events for the signal case we have not included decay channels with branching ratios smaller than 4%.
For this analysis, we have required the leptons to have at least p T > 25 GeV and |η| < 2.5, which removes the contribution of charged leptons produced in the forward and backward regions.We also required that the magnitude of the missing transverse energy to be at least 50 GeV.We then scanned over the threshold of selection on the transverse momentum of the charged lepton from 25 GeV to 1000 GeV and computed the significance for each selection threshold, including a 5% uncertainty on the background events.The results are shown in figure 9. We can see that the signal significance can easily reach 5 for the µ ± τ ∓ and µ ± µ ∓ channels at 10 TeV.The sensitivity reach for the µ ± e ∓ channel is not so promising except in BP2.

Detector level
We now turn into a realistic analysis of the ϕ + ϕ − channel at the detector level in the 2 muons plus missing energy final state.The basic kinematic observables are shown in figure 10 while the cutflow table for the event selection is displayed in Table VI.First, we require events to satisfy E miss T > 100 GeV.We also veto events that contain one isolated electron with p T > 10 GeV and |η| < 2.5, one tau lepton with p T > 25 GeV and |η| < 2.5 and one photon with p T > 100 GeV and |η| < 2.5.Since the two singly-charged scalars decay exclusively into muons, we require the presence of exactly two muons with p T > 30 GeV and |η| < 2.5.The two muons are required to have opposite electric charge.The ZZ backgrounds can be further suppressed by requiring that the invariant mass of the two muon system to be larger than 100 GeV.After this selection, about 50% of ZZ background events are removed while the contribution of W W/W Z backgrounds remains unchanged.The stransverse mass (M T 2 ) is a useful kinematic variable that can be employed to reduce the contributions of the W W backgrounds [97][98][99].Denoting the two charged lepton momenta by p T , and the missing particle momenta as q (1)

T and q
(2) T , the M T 2 variable is calculated as follows: T +q (2) T , q T ), M T (p T , q (2) where q (1) T and q (2) T are combined to form the total missing momentum.The calculation of the M T 2 also relies on the test mass for the invisible particle, which we choose to be equal to zero.For this choice, the M T 2 distribution will have an end point at around the mass of the parent particle, which can be seen clearly in the top right panel of figure 10 where the end of point of M T 2 is at around 1.25 TeV (BP1), 2.5 TeV (BP2), and 3.75 TeV (BP4).There are exceptions for the SM backgrounds since we have contributions of VBF diagrams that affect the expected behavior in s-channel mediated processes.We further require M T 2 > 100 GeV, which is enough to suppress the W W/W Z and ZZ backgrounds by 99% and 2%, respectively.To further reduce the contribution of ZZ backgrounds, we also impose the condition ∆ϕ(⃗ µ 1 , ⃗ p miss ) > 0.5 where ⃗ µ 1 is the 3-momentum of the leading muon.This selection reduces the yields of the two backgrounds by about 50%.Finally, we impose that M ℓℓ falls in the range ]200, 5000] GeV, which dramatically reduces the number of events for the ZZ backgrounds.The final selection efficiency for the signal events is of order 44%-56% (higher for higher scalar masses).
We finally calculate the signal significance in the signal region defined as the last step of the selection in Table VI.In this we assume that the number of backgrounds is the same for all the luminosities (n b ≡ N W W + N ZZ = 1.4), and we take into account two assumptions on the background uncertainty, i.e. assuming that δ = 0% (eq.26) and δ = 20% (eq.27).The results are shown in figure 11 where we can see that luminosities of about 200-400 fb −1 are enough to discover or exclude our signal benchmark points, therefore competing with or even outperforming the sensitivity reach of the HL-LHC.The case of the pair production of doubly-charged scalars is much more simpler since they decay predominantly into µµ with branching ratios of 93.6%-99.7%except in the case of BP5 where the decay branching ratio into muons is 63.1%.Using similar arguments regarding mass choices in our benchmark points, we find that BP1, BP2, and BP4 will lead to similar consequences since they correspond to m κ = 1.25 TeV.Therefore, in this channel, we will analyse the sensitivity reach in the following benchmark points: BP1, BP3, and BP5.The we also calculate the uncertainty on the number of events and the efficiency defined as ε = N i /N i−1 where N k refers to the number of events after the selection step k.All the events are normalized to 10 ab −1 of Luminosity.
number of signal and background events is defined as It is clear that we can easily reach signal significance larger than O(1) given the smallness of the associated background cross sections.We do not include the backgrounds from vectorboson fusion since they involve missing energy and extra leptons, and therefore, their inclusion would require a dedicated analysis strategy, which we will discuss in the next section.We first select events that contain exactly four muons with transverse momentum of at least 25 GeV and pseudorapidity smaller than 2.5.We also require that muons with opposite sign to be combined with each other to form a Z-boson candidate.
To further reduce the background contribution, we require that the invariant mass of these Z-boson candidates to be larger than 100 GeV (which kills most of the Z-boson contribution to the background).We then scan over the minimum transverse momentum of the muon for the range of 25-1000 GeV and compute the signal significance as defined in equations 26-27 including a 20% background uncertainty (dashed lines).The results are shown in figure 12, where we can see that the signal significance is very large for all the benchmark points.We conclude that the muon colliders will have a great potential to discover the Zee-Babu model, at least in channels involving doubly-charged scalars.We notice the importance of a dedicated analysis involving state-of-art Monte Carlo tools to comprehensively study the model at muon colliders, especially concerning reconstructions of scalars at muon colliders and their characterization.This will be discussed in the next section.

Detector level
We now turn to the analysis at the detector level using the same technical setup discussed in section VI A. The contribution of the ZZ production through VBF can be drastically reduced   by imposing cuts on missing energy (in case of charged-current production) and requiring exactly four muons (in case of neutral-current production).The cutflow table for this channel is shown in Table VII.We first select events if the total transverse missing energy is smaller than 100 GeV.Such a choice is motivated by the fact that the signal process does not involve any missing energy besides some minor contribution from misidentification of muon candidates.This selection step reduces the contribution of VBF ZZ by 93% while slightly affecting both the ZZ contribution through s-channel and the signal events for all the benchmark points.We then apply a number of vetoes on several objects like electrons with p T > 15 GeV and |η| < 2.5, τ leptons with p T > 25 GeV and |η| < 2.5, photons with p T > 100 GeV and jets with p T > 30 GeV and |η| < 2.5.We then require the existence of four charged muons forming two pairs of same electric charge.Those pairs will be used to form doubly-charged scalar candidates.The invariant mass of these candidates is required to be larger than 1000 GeV.The final efficiency for the backgrounds is small (0.09 for s-channel ZZ and 5.5 × 10 −4 for VBF ZZ).On the other hand, the signal efficiency is > 50%.To check the ability of our analysis in reconstructing the invariant mass of the κ candidates, we display the invariant mass of the system formed by same sign same flavor (SSSF) muons in figure 13 for BP1 (red), BP3 (blue) and BP5 (green).We can see from that the peak of the invariant mass distribution is centered around the mass of the doubly-charged scalar candidate, and the width of the distribution grows with the mass reflecting the full off-shell effects taken into account in our simulation.Given the large number of events that survive the final selection, we can conclude that the discovery reach for the case of κ ++ κ −− is even much more promising than the case of ϕ + ϕ − production.

VII. SUMMARY AND CONCLUSIONS
In this work, we have studied the phenomenology of the Zee-Babu model at future muon colliders.This model provides a natural explanation of the smallness of neutrino mass through two-loop radiative neutrino mass generation mechanism and with the exchange of two new SU (2) L singlet scalars: a singly-charged and a doubly-charged scalars.After studying the impact of the constraints from neutrino oscillation data, lepton flavour violation, and lepton flavour universality tests, we have selected five phenomenologically viable benchmark points which depend on the masses of the charged scalars of the model.We then analyzed the production rates for various channels that include charged leptons (both FC and FV channels), singly-charged scalars and doubly-charged scalars for the five benchmark points.A special attention was given to channels that not only probe the Yukawa-type couplings connected to the lepton sector but also channels that are sensitive to the nature of the scalar potential, in particular, the lepton-number violating coupling.These channels involve n ≥ 4 scalars and are enjoying almost a background-free environment even that the production cross sections are of order 10 −3 -10 0 fb.In particular, the production of four doubly-charged scalars lead to a golden channel consisting of eight highly-energetic muons.We finally analyzed the signal-tobackground for the pair production of two charged scalars: ϕ + ϕ − and κ ++ κ −− at √ s = 10 TeV both at the parton level and the detector level.For the pair production of singly-charged scalars (ϕ + ϕ − ), we have analyzed three channels, i.e. e ± µ ∓ +E miss T , µ ± τ ∓ +E miss T and µ ± µ ∓ +E miss T .For the pair production of doubly-charged scalars, we have analyzed the production of µ + µ − µ + µ − and we found even more promising sensitivity.

2 .
FIG. 2. Examples of Feynman diagrams for the pair production of charged leptons at muon colliders.The Drell-Yan type diagram on the left contributes only for the case of flavour-conserving case while the t-channel diagram on the right contributes to both the flavour conserving (α = β) and the flavour violating (α ̸ = β) production channels.
FIG. 3. Production cross section in femtobarn of FC lepton pair production as function of the centerof-mass energy for the muon channel (left) and electron channel (right).The results are shown for m κ = 1250 GeV (red), m κ = 2500 GeV (green) and m κ = 3750 GeV (blue).The results for m κ = 2500GeV and m κ = 3750 GeV for µµ production on the left panel are equal.We provide more details about this on the text.

FIG. 4 .
FIG.4.Production cross section in femtobarn of FV lepton pair production as function of the centerof-mass energy for the µ ± τ ∓ channel (left) and e ± µ ∓ channel (right).The results are shown for BP1

FIG. 9 .
FIG. 9.The signal significance as a function of the cut on the transverse momentum of the charged leptons for BP1 (left panel), BP2 (middle panel) and BP4 (right panel).Here we show the sensitivity for the µ ± e ∓ + E miss T

3 NeventsFIG. 10 .
FIG. 10.Kinematic distributions for the two muons plus missing energy final state at √ s = 10 TeV and L = 10 ab −1 .From top left to bottom right we show the transverse momentum of the leading muon (p T,µ 1 ), the transverse missing energy (E miss T ), the stransverse mass (M T 2 ), the invariant mass of dimuon system (M ℓℓ ) and the difference in the azimuthal angle between the leading muon and the missing energy momentum (∆ϕ(µ 1 , p miss )).The backgrounds are stacked on the top of each other where we show ZZ (yellow) and W W/W Z (cyan).The signal is shown for BP1 (red), BP2 (green), and BP4 (blue).

FIG. 11 . 3 √ 4 BP3µ 4 BP5FIG. 12 .
FIG. 11.Signal significance for the discovery reach of the singly-charged scalar in the µ + µ − + E miss T final state as a function of the luminosity.Here we show the results for BP1 (red), BP2 (blue), and BP4 (green).The results shown assuming 0% and 20% uncertainty on the background yields in dashed and solid lines, respectively.

TABLE III .
Constraints from loop-level lepton flavour violating interactions and anomalous magnetic moments.Experimental limits are taken from

TABLE IV .
The benchmark points used in this study.

TABLE V .
Predictions of the model of the neutrino masses and mixings, and lepton-flavor violating decay branching ratios.The neutrinoless double beta decay parameter is defined as m ββ ≡ i U 2 ei m i .Predictions of the cLFV modes that are just below the current bound and will be fully tested in the next upgrades are highlighted with boldface.

TABLE VII .
Same as inTable VI but for the case of the four muon final state.