Full parameter scan of the Zee model: exploring Higgs lepton flavor violation

We study the general Zee model, which includes an extra Higgs scalar doublet and a new singly-charged scalar singlet. Neutrino masses are generated at one-loop level, and in order to describe leptonic mixing, both the Standard Model and the extra Higgs scalar doublets need to couple to leptons (in a type-III two-Higgs doublet model), which necessarily generates large lepton flavor violating signals, also in Higgs decays. Imposing all relevant phenomenological constraints and performing a full numerical scan of the parameter space, we find that both normal and inverted neutrino mass orderings can be fitted, although the latter is disfavored with respect to the former. In fact, inverted ordering can only be accommodated if $\theta_{23}$ turns out to be in the first octant. A branching ratio for $h \to \tau \mu$ of up to $10^{-2}$ is allowed, but it could be as low as $10^{-6}$. In addition, if future expected sensitivities of $\tau\to \mu\gamma$ are achieved, normal ordering can be almost completely tested. Also, $\mu e$ conversion is expected to probe large parts of the parameter space, excluding completely inverted ordering if no signal is observed. Furthermore, non-standard neutrino interactions are found to be smaller than $10^{-6}$, which is well below future experimental sensitivity. Finally, the results of our scan indicate that the masses of the additional scalars have to be below $2.5$ TeV, and typically they are lower than that and therefore within the reach of the LHC and future colliders.


Contents
1 Introduction

The scalar potential
The following analysis is similar to the ones performed for 2HDMs, see e.g. refs. [72][73][74][75]. One can start in a generic basis, where both Higgs scalar doublets Φ 1 and Φ 2 take VEVs denoted by v 1 and v 2 , respectively. Then, one can rotate to the Higgs basis, where only H 1 takes a VEV denoted as usual as v = v 2 1 + v 2 2 246 GeV. The rotation is given by the following transformation [73]: where tan β ≡ v 2 /v 1 and the short-hand notations s x ≡ sin x and c x ≡ cos x. 3 We will also use t x ≡ tan x. In the Higgs basis, the doublets take the form: where the CP-even mixing is given by that needs to be sufficiently close to zero (i.e. the alignment limit) to give rise to a SM-like Higgs boson [72].

The lepton sector
As we will see, in order to describe leptonic mixing, both Higgs scalar doublets must couple to the charged leptons, and thus, we are considering a type-III 2HDM, see e.g. ref. [51]. The most general Yukawa Lagrangian in the generic basis, where both Higgs fields take VEVs, reads where L = (ν L , e L ) T and e R are the SU(2) lepton doublets and singlets, respectively, and L ≡ iσ 2 L c = iσ 2 CL T with σ 2 being the second Pauli matrix. Due to Fermi statistics, f is an antisymmetric Yukawa matrix in flavor space (i.e. f αβ = −f βα ), while Y 1 and Y 2 are completely general complex Yukawa matrices. Furthermore, the charged-lepton masses are given by (2.16) Note that we will work in the basis where m E is diagonal with real and positive elements m e , m µ , and m τ . Moreover, Y 2 will be a general complex matrix and Y 1 can be expressed completely in terms of m E and Y 2 using eq. (2.16).
In the Higgs basis, we can rewrite eq. (2.15) using eq. (2.1) as H 2 e R +Lf Lh + + H.c. (2.17) Without loss of generality, rotating the lepton doublets L α and the lepton singlets e α R by the same phase (so that m E remains diagonal and positive), three phases from f can be removed. However, note that the phases from Y 2 cannot be removed by lepton field redefinitions.
In the mass basis (also for massive neutrinos), using eq. (2.17), the most general leptonic Lagrangian reads 2c β e R A + H.c. (2.18) We define the following effective couplings for the neutral Higgs fields h 0 = (h, H, A), which will turn out to be useful for CLFV processes: (2.20) and for the charged Higgs fields h c = (h + 1 , h + 2 ): One can observe that g 1 h 0 is flavor conserving and proportional to m E /v.
Using the three definite neutrino masses m 1 , m 2 , and m 3 , we also define the two linearlyindependent neutrino mass squared-differences ∆m 2 21 ≡ m 2 2 − m 2 1 and ∆m 2 31 ≡ m 2 3 − m 2 1 , known as the small and large mass squared-differences, respectively, where the sign of ∆m 2 31 is still unknown. The case ∆m 2 31 > 0 is generally referred to as 'normal neutrino mass ordering' (NO), whereas the case ∆m 2 31 < 0 is known as the 'inverted neutrino mass ordering' (IO). Note that using neutrino oscillation experiments, it is not possible to determine φ 1 and φ 2 nor the absolute neutrino mass scale. The most up-to-date bestfit values from global analyses of the ordinary neutrino oscillation parameters (i.e. the leptonic mixing parameters and the neutrino mass-squared differences) are θ 12 34 • , θ 13 8.5 • , θ 23 42 • for NO and θ 23 50 • for IO, δ 1.5π, ∆m 2 21 7.5 · 10 −5 eV 2 , and ∆m 2 31 2.5 · 10 −3 eV 2 for NO and ∆m 2 31 −2.5 · 10 −3 eV 2 for IO [59][60][61]. Similar values can also be found in refs. [77,78]. Note that both the first and second octants of θ 23 are allowed from the global analyses [60,61] for both orderings with a mild preference of the first (second) octant for NO (IO).
In addition to the ordinary neutrino oscillation parameters, the following effective neutrino parameters appear naturally in different contexts [76] where m ee (or m 0ν2β ) is the effective electron neutrino mass parameter that could be measured in neutrinoless double beta decay (0ν2β) experiments [79,80] (see also ref. [81] for a recent review), m νe (or m β ) is the effective neutrino mass parameter measured in (single) beta decay experiments [82], and finally, m i is the sum of the three neutrino masses, which, in the future, could be determined by cosmology, but at present it is only restricted by an upper bound, see e.g. refs. [83,84].

Neutrino masses in the Zee model
As can be seen from the potential and the Yukawa Lagrangian of the Zee model, eqs. (2.3) and (2.18), respectively, in order to have LNV and therefore neutrino masses, we need the simultaneous presence of Y 1 , Y 2 , f , and µ. In the Zee model, the one-loop diagram shown in figure 1, where the charged scalars h + 1 and h + 2 run in the loop, generates neutrino masses. The complete neutrino mass matrix is then given by (see e.g. ref. [54]) with ϕ being the mixing angle for the charged scalars given in eq. (2.8). Therefore, in the Zee model, due to loop and chiral suppression, the new physics scale can be light. Assuming f eµ = 0, neglecting m e m µ , m τ , and keeping only the term proportional to m µ in the 3-3 element, 4 we obtain the following (symmetric) Majorana mass matrix (2.32) Note that in a simpler scenario, where one neglects terms proportional to m µ , one neutrino will be massless. On the other hand, taking terms proportional to m µ into account, all neutrinos will obtain masses.
In our analysis, we assume zero Yukawa couplings in the e-µ sector, i.e. we assume f eµ = 0 and Y µµ We assume all these Yukawa couplings to be complex except for f eτ , f µτ , and Y eτ 2 (which does not enter in neutrino masses), see the discussion in section 2.2. Thus, the Yukawa couplings will constitute eleven free real parameters in the numerical scan that will be described in section 4.1.
In order to obtain correct mixing angles, we need both Y τ µ 2 and Y τ e 2 different from zero, as they enter in the 1-2 submatrix of eq. (2.32). Therefore, it is clear that reproducing the leptonic mixing angles correctly will imply restrictions on Br(h → τ µ), Br(h → τ e), and other LFV processes. In fact, from this argumentation, a lower bound on the product Br(h → τ µ) · Br(h → τ e) (in addition to an upper bound from other CLFV processes) is expected.

The (minimal) quark sector
Although the Zee model only deals with the lepton sector, the SM Higgs scalar doublet needs to couple to the SM quarks, like tops and bottoms, in order to be observed via its production and decay modes at the LHC [85]. In the generic basis, the most general Lagrangian in the quark sector is given by where Q = (u L , d L ) T are the SU(2) quark doublets, u R and d R are the SU(2) quark singlets, andΦ i ≡ iσ 2 Φ * i (i = 1, 2). However, flavor violation in the quark sector is severely constrained (see e.g. ref. [16]), so we will assume the simplest scenario in which Y d2 = Y u2 = 0. Then, we can use the basis, where the up-type quark mass matrix is diagonal. Furthermore, we assume the Yukawa couplings Y d1 and Y u1 to be Hermitian. The masses for the third generation quarks are given by Therefore, the interactions of the physical neutral Higgs fields with quarks are given by where the corresponding Feynman rule for the CP-odd scalar A includes a γ 5 . Note that if we had taken the couplings to quarks as general as possible, including the first generation, there would have been other phenomenological implications. In particular, related to neutrino masses, there would have been be new contributions to neutrinoless double beta decay and new universality and non-standard neutrino interactions with matter, stemming from interactions of the charged scalars h + 1 and h + 2 . However, when naturality constraints are imposed on the Yukawa couplings to the leptons and the up and down quarks, see eq. (3.4), these contributions are subdominant. We will therefore only discuss the universality constraints and the non-standard neutrino interactions generated through leptonic interactions, see section 3.4, and we will only consider the contributions to neutrinoless double beta decay mediated by W bosons, i.e. the contributions from the light neutrinos.

Stability of the potential
A Hamiltonian in quantum mechanics has to be bounded from below, which requires the quartic part of the scalar potential in eq. (2.3) to be positive for all values of the fields and for all scales. Then, if two of the three fields H 1 , H 2 , and h vanish, one immediately finds For a general 2HDM potential with λ 6 = λ 7 = 0, it has been shown in ref. [86] that the additional necessary conditions are However, when λ 6 , λ 7 = 0, it has been shown that in addition to the previous conditions, the following condition [87] 2 |λ 6 + λ 7 | < λ 1 + λ 2 2 is both necessary and sufficient to ensure stability of the potential. Other stability conditions for similar potentials are discussed in refs. [45,88,89]. In this work, due to the large number of parameters, we will set λ 4 = λ 7 = λ 8 = λ 9 = λ 10 = λ h = 0, since they do not significantly impact phenomenology, even though their presence is expected to somewhat open the allowed parameter space. Thus, the four free Higgs couplings are λ 1 , λ 2 , λ 3 , and λ 5 , which we will treat as free real parameters, while λ 6 is a derived parameter that can be computed from eq. (2.13). In the numerical scan (see section 4.1), we impose the conditions from eqs.

Naturality and perturbativity
There are naturality and perturbativity constraints on the Yukawa couplings and on the quartic and trilinear couplings of the potential. In order not to have large fine-tuned cancellations between the different Yukawa couplings (see e.g. refs. [20,23] One can also obtain an upper bound on µ, which contributes to the scalar masses. In fact, using eq. (2.9), it is clear that naturality demands that µ √ we can also derive a natural upper bound using the 125 GeV Higgs boson [7,8], due to the fact that µ contributes at one-loop level to its mass. The relevant coupling of the light Higgs boson to the charged scalars induced by µ is We demand that the one-loop contribution to the Higgs mass fulfills δm h /m h κ, where we choose κ = 1 (10), which corresponds to no (10 %) fine-tuning. 5 Neglecting logarithms and factors of two in the Higgs self-energies, we obtain Taking s β−α ∼ 1, we find an upper bound of 1.5 (15) TeV for κ = 1 (10). In addition, we impose that all the quartic couplings are perturbative: 3.3 Charged lepton flavor violation and electric and magnetic moments

Trilepton decays
The presence of the second Higgs doublet gives rise to tree-level trilepton decays i → j k l . 6 The ratio of branching ratio reads The fine-tunings in the Higgs mass squared, which is the relevant parameter in the Lagrangian, would be 1 % (100 %) for κ = 1 (10). 6 At one loop and two loops, there are dipole contributions which dominate the rate. However, these are strongly bounded by τ → µγ and µ → eγ, see section 3.3.2. Also, box diagrams are very suppressed, see ref. [90].
where G F = g 2 /( √ 2v 2 ) 1.166·10 −5 GeV −2 is the Fermi coupling constant and the Wilson coefficients D PP (P, P = L, R) are given by the coherent sum of the contributions from the neutral Higgs fields. In addition, ξ = 1/2 (1) when there are two (no) indistinguishable particles in the final state. We are interested in tau decays. 7 In this case, for τ → µµμ, the Wilson coefficient D LL is given by where h 0 = (h, H, A). Similarly, for τ → eµμ, one can simply substitute Y µτ 2 ) * and conjugating the vertex factors.
As expected, these tree-level processes do not restrict the parameter space as much as τ → µγ, τ → eγ, or µ → eγ does, since they always involve a muon mass suppression (squared) and are therefore irrelevant. In table 5, the upper bounds used in the numerical scan for the various observables are presented.

3.3.2
i → j γ decays One of the most constrained CLFV process is the radiative process i → j γ with i being the physical charged leptons e, µ, and τ . This process always arises at loop level and it can be viewed as stemming from an effective operator of the form [21] where Λ is the scale of new physics, σ µν = i[γ µ , γ ν ]/2, and F µν = ∂ µ A ν − ∂ ν A µ is the electromagnetic field strength tensor with A µ being the photon field. It is useful to define 11) and similarly, C L = C † R . Then, it follows that We use the expressions for the Wilson coefficients given in refs. [13,15,16], adapted to the Zee model. We also include the two-loop Barr-Zee contributions as given in ref. [91]. At one-loop level, the dominant contribution for τ → µγ reads Other processes, like µ → eee, are absent at tree level, since we assume Y µe 2 = Y eµ 2 = 0. Also, tree level contributions to τ → µeē and τ → eeē are suppressed by me. and similarly, we obtain C 0 R with the replacement Y αβ 2 → (Y βα 2 ) * . For the charged scalars, we find (using i |U τ i | 2 ≈ 1) where the contribution to C + R is zero, since we assume f eµ = 0. The total contribution to the Wilson coefficient is and similarly for C R . For µ → eγ, the dominant contributions, proportional to m τ , are given by the neutral Higgs fields and τ running in the loop Note that we also add the contributions from f to µ → eγ, which are proportional to m µ , with ν τ running in the loop: The contribution in eq. (3.19) strongly constraints the antisymmetric Yukawa coupling f of the singly-charged scalar singlet.

Electron and muon electric dipole and anomalous magnetic moments
When flavor is conserved, anomalous magnetic moments (AMMs) are generated. The dominant contributions are given by loops with neutral Higgs fields and tau leptons, similar to those in eqs. (3.16) and (3.17). It can be defined as [15] Then, using eq. (3.11), the muon AMM is given by Similarly, the electron AMM is obtained by replacing the indices µ → e everywhere. For the electron AMM, there is no disagreement between theory and experiment, and in fact, it represents one of the most precisely measured quantities in all of physics with an experimental 90 % C.L. upper bound on possible new physics contributions of 2 · 10 −12 [76]. On the other hand, for the muon AMM, there is an experimental deviation from the theoretical prediction of the SM, i.e. |δa µ | = (2.88±0.80)·10 −9 [76]. Unfortunately, neither a 2HDM nor the Zee model can accommodate this discrepancy [23].
If there is CP violation, electric dipole moments (EDMs) are generated, which can be defined as The experimental 90 % C.L. upper bounds for the muon and electron EDMs are 1 · 10 −19 e cm [76] and 8.7 · 10 −29 e cm [92], respectively (see also table 5). Therefore, for the muon EDM, this bound is almost non-existing, whereas the electron one strongly constrains the imaginary parts of Y τ e 2 and Y eτ 2 (in our scenario, Y eτ 2 is assumed to be real for simplicity).

µe conversion in nuclei
Finally, µe conversion is a very interesting process because the sensitivity in the next generation experiments is expected to increase by around four orders of magnitude. As the dipole contribution is already heavily constrained by µ → eγ, we can use the current limits to restrict the monopole photonic contribution. In the limit, where the transferred momentum is zero, the conversion rate, relative to the muon capture rate, can be expressed in gold as [93] Cr where V p Au = 0.0974 m 5/2 µ , Γ Au = 8.6 · 10 −18 GeV [93], and the vector current coefficient is given by [23] Similarly, g RV is obtained by replacing Y αβ 2 → (Y βα 2 ) * in eq. (3.25). The currently best experimental limit is Cr(µ → e) Au < 7·10 −13 at 90 % C.L. by the SINDRUM II experiment at PSI [94].
For titanium, the best experimental limit is Cr(µ → e) Ti < 4.3·10 −12 at 90 % C.L. [94]. In the future, PRISM/PRIME [68,69] expects to achieve a sensitivity of O(10 −18 ). For aluminium, the future experimental sensitivity is Cr(µ → e) Al < 6 · 10 −17 at 90 % C.L. in the Mu2e experiment at Fermilab [64][65][66] and of O(10 −17 ) in the COMET experiment at J-PARC [70,71]. Mu2e may also run at Project X using either Al or Ti [67] with an expected sensitivity of O(10 −19 ). In the numerical analysis, we will discuss how the Zee model will be constrained if the expected sensitivity of this process is achieved in the future.

Leptonic interactions of the charged scalars
Now, we study four-lepton interactions of the singly-charged scalar mass eigenstates h + 1 and h + 2 that involve two charged leptons and two neutrinos. 8 These give rise to muon and tau decays into lighter charged leptons and neutrinos as well as to non-standard neutrino interactions (NSIs), see e.g. refs. [95,96].
At tree level, we can integrate out h + 1 and h + 2 , see e.g. ref. [95]. Therefore, using eq. (2.18) and expanding to first order in p 2 /m 2 where we have defined the effective Yukawa couplings and the effective masses (3.28) Note that for the definition ofM 2 12 we have used the charged-scalars mixing s 2ϕ as defined in eq. (2.8).

Universality
The second operator in eq. (3.26), which is second order in f , i.e. |Y eff 2 | 2 ∝ f † f , couples to the left-handed leptons, like charged currents in the SM. This implies that it interferes constructively with the W boson. In the SM, the Fermi constant extracted from muon decay G SM µ and the one extracted from hadronic decays G SM β are tested to be equal with great precision. The presence of the charged scalars modifies the muon decay rate [97,98], which implies that G SM β = G SM µ = G Zee µ , where G Zee µ is the Fermi constant from muon decay in the Zee model. Therefore, we have whereM 2 is defined in eq. (3.28).
In the SM, unitarity of the quark mixing matrix V holds to great precision. In our scenario, as we assume f eµ = 0, we also have that V is unitary up On the hand, other leptonic decays may not be universal (in the SM, they are mediated by gauge interactions and are therefore universal). The ratio of flavor violating decays can be tested among the different generations via the effective couplings given by (3.33) In our scenario the expressions (3.31)-(3.33) will generally, but not necessarily, be different from one, i.e. they deviate from the SM prediction.

Non-standard neutrino interactions
Apart from standard neutrino interactions (including neutrino oscillations), the new singlycharged scalar fields h + 1 and h + 2 introduced in the Zee model will induce NSIs at tree level. These NSIs are new LFV processes that are not allowed in the SM, but could be probed in future neutrino oscillation experiments, and are usually treated using an effective fourfermion operator.
Interestingly, the operators in the second line of eq. (3.26) violate lepton number. 9 Indeed, they involve the same combination of four leptons that appears inside the neutrino mass diagram, see figure 1, and thus, their coefficients are proportional to the same leptonnumber combination appearing in the neutrino mass formula, see eq. (2.32). Hence, they are subject to constraints from neutrino masses and therefore suppressed.
The operators in the first line of eq. (3.26) do not violate lepton number and, in principle, they give rise to NSIs that are not suppressed by neutrino masses. Applying Fierz identities one can express them in various ways. Using ref. [100], they can be written in the flavor basis with the usual NSI language as [101,102] where χ ρσ αβ and ε ρσ αβ are the canonical NSI parameters given by whereM 1 andM 2 are defined in eq. (3.28). For neutrinos propagating in ordinary matter, these operators induce the following matter NSI parameters (3.36) 9 Their gauge invariant EFT operators are, of course, dimension 7 [99], LLLeRΦ1,2.
In our scenario, the only relevant matter NSI parameters are χ m τ τ and ε m τ τ (cf. ref. [102]), since we assume Y µe 2 = Y eµ 2 = 0 and f eµ = 0. Now, we can derive an upper bound on χ m τ τ applying CLFV and HLFV limits. For illustration, let us impose the limits from Br(h → τ e) from table 1. Using a similar equation to eq. (3.48), but for the τ e channel, which depends on Y eτ 2 , we obtain which is below present and most probably future experimental sensitivity. Reproducing small neutrino masses and fulfilling other stronger CLFV constraints, e.g. constraints on µ → eγ and τ → µγ, imply that χ m τ τ and ε m τ τ are much smaller than 10 −8 , which is beyond any future experimental sensitivity. This can be seen in our numerical scan. Other limits and future prospects on NSIs can be found in refs. [103,104].
Finally, at a future neutrino factory, there could be source NSIs in the process µ → eν β ν α . In our scenario, the only relevant source NSI parameters χ µe αβ and ε µe αβ are those with tau neutrinos, i.e. χ µe τ τ and ε µe τ τ . However, these are also very small, at least below 10 −6 .

Higgs signals
In the Zee model, the couplings to SM particles are modified with respect to their SM values. For instance, the couplings to gauge bosons are and similarly, for g hZZ , g HZZ , and g AZZ , changing m 2 W → m 2 Z . Clearly, close to the decoupling limit, s β−α → 1, the light Higgs interactions are sufficiently SM-like [72]. As we will see in section 3.5.2, in order to have HLFV, we cannot be exactly at the decoupling limit, but it needs to be close enough to fulfill the bounds on the Higgs decays measured at the LHC [85]. The other Higgs couplings in the Zee model are modified as in eq. (2.18) (see also eq. (3.48)) for leptons, eq. (2.34) for quarks, and eq. (3.42) for photons.
The Higgs results at the LHC are usually given in terms of the global signal strength defined as where σ X (h) is the cross section for the production mode X and Br(h → Y ) is the Higgs branching ratio for the decay mode Y . By definition, in the SM, µ SM if = 1 for all production modes i and decay channels f . At the LHC, there are four production modes available for the Higgs boson, where the dominant one is gluon-gluon fusion (ggF), mainly through a top loop. The subdominant ones are vector boson fusion (VBF), associated production with a vector boson V h (where V = W, Z), and the associated production with a top-quark pair tth. The production modes are usually grouped into two effective modes according to ggF + tth and VBF + V h. We consider the five decay channels, where a signal has been detected, namely γγ, W W * , ZZ * , bb, and ττ . For instance, for gg → h → bb, the signal strength is given by For the contribution to the χ 2 function (to be discussed in section 4) from the Higgs decay channels, we need to take into account correlations between different production modes. Thus, for each of the decay modes f = γγ, W W * , ZZ * , bb, ττ , the contribution to the χ 2 function is defined as

h → γγ decays
In the Zee model, the decay of the Higgs boson to two photons is modified by two factors. First, the couplings to gauge bosons and top quarks are changed, since we have two Higgs doublets. Second, there are new extra charged scalars couplings to the Higgs boson. In ref. [45], a study of h → γγ in the Zee model has been performed. However, λ 7 and λ 10 were set to zero. Thus, in the following, we will analyze this decay in our scenario.
The value of the h → γγ decay width in the Zee model with respect to the SM one is given by [109][110][111] where λ SH is the coupling of a charged scalar S with mass m S to the Higgs field, which is coming from a term in the potential of the form (H † H)(S † S). Note that we have used the modified couplings to tops and W W given in eqs. (2.34) and (3.38), respectively. Here, τ i ≡ 4m 2 i /m 2 H and A i (x) (i = 0, 1/2, 1) are loop functions: We need to compute the couplings to charged scalars λ SH in the Zee model. Since we are in the Higgs basis, the terms of the potential in eq. (2.3), which are coupling the Higgs boson to the charged scalars, are those involving λ 3 , λ 7 , λ 8 , λ 10 , and µ. Using the rotations to the mass basis, i.e. eqs. (2.7) and (2.12), the relevant interactions hh + 1 h − 1 and hh In the case when there is no mixing µ (ϕ) → 0, we obtain where the first equation agrees with eq. (F1) in ref. [72]. Note that, in our scenario, only the terms proportional to λ 3 and µ in eq. (3.46) will contribute, as we set the other couplings to zero.
Note that for a type-III 2HDM there is an upper bound on Br(h → τ µ) · Br(h → τ e) from combining the rates of µ → eγ and µe conversion [23] (which currently saturates the bound), as all combinations of couplings relevant to these HLFV processes enter in CLFV with tau leptons running in the loop, see secs. 3.3.2 and 3.3.4: (3.51) In the Zee model, as we will see, reproducing the leptonic mixing angles correctly implies that there are also lower bounds on the HLFV processes.

Scan of the parameter space
In order to study the large parameter space of our scenario of the Zee model and to be able to investigate how large CLFV and HLFV processes can be, we perform a full numerical scan using the software MultiNest [113][114][115]. MultiNest is a Bayesian inference tool that uses so-called nested sampling and especially suitable when there are possibly several maxima in the parameter space. It is designed to determine the Bayesian evidence, but as a byproduct, it also yields the posterior distribution that is relevant for a Bayesian analysis. Nevertheless, it also maximizes the likelihood, which is relevant for a frequentist analysis. We are interested in the maximization of the likelihood and we perform a fully frequentist analysis. We scan over all free parameters in the model, in total 19 real parameters, which are given in table 2 together with their chosen allowed parameter ranges. The plots, including best-fit points and 1σ and 2σ confidence regions, are produced using the graphical interface Superplot [116]. We impose the stability conditions on the scalar Table 2: Priors on the 19 free real parameters used in the scan. For Y τ τ 2 , Y τ µ 2 , Y τ e 2 , and Y µτ 2 , we scan real and imaginary parts independently, while f µτ , f eτ , and Y eτ 2 can be assumed to be real without loss of generality. We use logarithmic priors for all parameters except for tan β, where uniform priors are used. The quantity that is maximized is the likelihood L, which is equivalent to minimizing the χ 2 function: χ 2 = −2 ln L. We assume Gaussian likelihoods, and thus, the contribution to the χ 2 function from the N signals is given by is the theoretical prediction (experimental measurement) of the observable i with the respective standard deviation σ i . For s 2 23 and δ, we use their combined two-dimensional distribution function (therefore including their correlation) given in refs. [60,61]. For the M upper bounds, we use where B j is the experimental upper limit at 1σ significance level of the observable j.
According to the definitions for the standard Gaussian distributions assumed for the observables, the 90 % C.L. and 95 % C.L. upper limits are normalized (i.e. divided) by the factors 1.645 and 1.949, respectively. The limits for the signals are rescaled in the same way. The choice of eq. (4.2) for the χ 2 function is made so that we do not penalize deviation from zero when this is not supported by data.
We take into account all the relevant bounds and observables described in section 3, except for the muon AMM discrepancy with respect to the SM, which cannot be accommodated in our present scenario of the Zee model. Thus, adding eqs. (4.1) and (4.2), the total χ 2 function reads χ 2 = χ 2 signals + χ 2 bounds .
Thus, in our scenario of the Zee model, neutrino masses and leptonic mixing can be accommodated in both NO and IO. However, IO is disfavored compared to NO. We also note that if θ 23 will turn out to be in the second octant, IO cannot be accommodated.
In figure 2, we present the contributions to χ 2 min from the observables in both NO and IO. Note that we do not show the contributions from the upper bounds, since these Observable Upper bound 4.4 · 10 −8 [76] Br(τ − → e − γ) 3.3 · 10 −8 [76] Br(µ + → e + γ) 4.2 · 10 −13 [62] Br(τ − → µ − µ + µ − ) 2.1 · 10 −8 [76] Br(τ − → µ − µ + e − ) 2.7 · 10 −8 [76] Cr(µ → e) Au 7 · 10 −13 [94] |δa e | 2 · 10 −12 [76] |d µ /e| [cm] 1 · 10 −19 [76] |d e /e| [cm] 8.7 · 10 −29 [92]   are always satisfied, and thus, the corresponding contribution to χ 2 min is exactly zero. Note that we present the combined contribution from s 2 23 and δ to χ 2 min . In IO, the largest contributions stem from s 2 12 and s 2 23 +δ, although they are all within 3σ of their experimental values. In NO, the corresponding contributions are small. These observables account for the fact that the fit in IO is much worse than in NO. All other observables in IO are of the same size as those in NO and within 2σ. In addition, we perform a run assuming Y µτ 2 = 0, which renders one neutrino massless, see eq. (2.32). In this case, we find that IO is significantly better than NO, which is in agreement with the results of ref. [54], and the value of χ 2 min in IO is of the same size as for Y µτ 2 = 0 (i.e. quite large), whereas in NO, the difference is more than one order of magnitude with χ 2 min ∼ O (100). Thus, in this scenario, both NO and IO are basically excluded. Moreover, there are non-negligible contributions to χ 2 min from the Higgs signals, especially h → ZZ. Regarding h → γγ, we find that it is around the SM value. Note that in this case we choose some of the couplings of the scalar potential to be zero, see eqs. (3.46), so its value can be modified by turning them on. The values for the universality parameters g ij , see eqs. (3.31)- (3.33), are always very close to one, even closer than the experimental values, and compatible with observations at 2σ. However, these also give non-negligible contributions to χ 2 min . In the following, we present the allowed regions for the most interesting parameters and observables. In figure 3, we plot the leptonic mixing parameters s 2 12 and s 2 23 for NO (left panel) and IO (right panel), where one can clearly see that the fit is very good for NO, while for IO it crucially depends on the octant of θ 23 . In fact, the fitted value of s 2 12 would be 5σ away from its experimental best-fit value if θ 23 lies in the first octant. For NO, it is clear that both octants are viable.
In figure 4, we plot the allowed regions of the branching ratios Br(h → τ µ) and Br(τ → µγ). For µ = 0 (left panel), only upper bounds on the CLFV and HLFV processes exist and these are compatible with the results of refs. [22,23]. For µ = 0, one can observe that there are lower bounds on the CLFV and HLFV processes in both NO (middle panel) and IO (right panel). That is, reproducing neutrino masses implies that CLFV and HLFV on Br(τ → µγ) of O(10 −9 ) [63], which would significantly probe a substantial part of the allowed parameter space of IO and almost the complete allowed region of NO. This is one of the most interesting results of our work. Similarly, if sensitivities of Br(h → τ µ) at future colliders reach 10 −4 (10 −5 ), NO (IO) will be tested at 68 % C.L.
In figure 5, we show the allowed regions of the branching ratios Br(h → τ e) and Br(h → τ µ) for NO (left panel) and IO (right panel). No correlation exists for µ = 0, while there is a strong correlation for µ = 0, stronger for IO than for NO. In general, we find that Br(h → τ e) 10 −2 Br(h → τ µ) or even lower for both orderings. Therefore, observations of h → τ e will be considerably more challenging than for h → τ µ. In figure 6, we display the allowed regions of the µe conversion rate in gold and the branching ratio Br(h → τ µ) for NO (left panel) and IO (right panel). As discussed in section 3.3.4, next generation experiments, using aluminium and titanium, are expected to achieve an improved sensitivity of up to about four orders of magnitude, maybe reaching O(10 −19 ) [66]. The µe conversion rates of these materials are of the same order of magnitude as that of gold, 11 and therefore, if a negative result is obtained, IO would be excluded, while there would still be a considerable allowed region for NO. Thus, there is a complementarity between Br(τ → µγ) and µe conversion. On the other hand, the sensitivity of Br(τ → eγ) is expected to reach about 3 · 10 −9 [63] and Br(µ → eγ) is expected to be improved by one order of magnitude, but these are not able to test the model as thoroughly as Br(τ → µγ) and µe conversion.
In figure 7, we present the allowed regions for the absolute values of the Yukawa couplings Y τ e 2 and Y τ µ 2 in both NO (left panel) and IO (right panel). We find that the regions are quite well defined, although they are larger in IO than in NO, and the allowed values are larger for Y τ µ 2 than for Y τ e 2 . Note that the value of Y τ µ 2 , which is a very relevant parameter since it controls the decays of the scalars into τ µ (together with Y µτ 2 ), is always larger than around 10 −3 (2.5 · 10 −3 ) in NO (IO).
The scale of neutrino masses is controlled by the charged scalar mixing angle ϕ, which is proportional to the trilinear coupling µ, see eq. (2.8), and the Yukawa couplings f and Y 2 . Therefore, in figure 8, we plot the allowed regions of |f µτ | and s 2ϕ for both NO (left panel) and IO (right panel). The range of |f µτ | is similar in both orderings. One can clearly see that s 2ϕ is close to zero for |f µτ | 10 −5 , while it grows very fast for |f µτ | 10 −5 , reaching values of 0.7 (0.4) for NO (IO). As expected, the naturality condition of eq. (3.6), which is  added to the χ 2 function, restricts µ to be smaller than about 3 (30) TeV for κ = 1 (10) at 2σ.
The Yukawa couplings f µτ and f eτ are always in the range [10 −7 , 0.1] and have similar allowed regions in both NO and IO. In fact, they are highly correlated with f µτ f eτ (0.1f eτ ) in NO (IO). Their allowed 1σ C.L. regions lie roughly below 10 −5 , suppressing all interactions mediated by the antisymmetric Yukawa coupling f of the singlycharged scalar singlet. Therefore, in order to describe neutrino masses and leptonic mixing, Y 2 should be much larger than f .
In figure 9, we show the allowed regions of the Higgs scalar mass differences m H − m A and m H − m h + 1 . These affect the size of the parameter T , see appendix A. We find similar results to those in ref. [125]. At 1σ C.L. and at low scalar masses, m H −m h + 1 is roughly equal to m H − m A , thus canceling the contributions to T . At 2 σ C.L., for the two cases µ = 0 , the smaller the s 2ϕ , see eq. (2.8), and therefore, the smaller the neutrino masses. 12 From eq. (3.49), we know that Br(h → τ µ) is proportional to 1/m 4 H , i.e., it decouples with the CP-even scalar mass. In figure 10, we display the allowed regions of the mass m H and the branching ratio Br(h → τ µ). For NO (left panel), m H 0.9 (1.7) TeV at 1 (2) σ C.L., while for IO (right panel), m H 0.7 (1.1) TeV at 1 (2) σ C.L. Therefore, if an extra CP-even scalar (and close by CP-odd and charged scalars, see figure 9) is observed below 0.9 (0.7) TeV, then one expects Br(h → τ µ) 10 −4 (10 −5 ) for NO (IO) at 1σ C.L. The heavy CP-even scalar H could also have sizable decays into τ µ, depending on the scalar spectrum. In figure 11, we show tan β as a function of s α for µ = 0 (left panel), NO (middle panel), and IO (right panel). Having Higgs boson decays close to the observed ones implies being close to the decoupling limit, i.e. s β−α → 1, and therefore, tan β and s α are strongly correlated. For massless neutrinos, tan β can reach values up to 15 and s α can approach zero. However, if neutrino masses are introduced, the value of tan β is severely constrained to smaller values in both orderings and the allowed range of s α is also reduced. The upper bound on tan β is clearly more severe for NO, where tan β is driven to the smallest possible values (i.e. below 0.5), while for IO, it can reach values up to 1.4. We have also performed a run forcing the value of the unphysical parameter tan β to be large, i.e. 40 tan β 50, and we find that the fit becomes significantly worse, leading to a value for the minimum of the χ 2 function of O(1000), and thus excluding this scenario.
In figure 12, we present the allowed regions of the effective mass parameter m ee that appears in neutrinoless double beta decay and the smallest neutrino mass (m 1 for NO and m 3 for IO) for NO (left panel) and IO (right panel). We can see that the smallest neutrino mass is several orders of magnitude smaller in IO than in NO. In fact, in IO, it can be massless, whereas this is not the case in NO. This is consistent with the fact that the fit in NO is bad when Y µτ 2 = 0, rendering one neutrino massless. In NO, we obtain |m ee | (4 − 5) meV, while in IO, |m ee | is one order of magnitude larger, i.e. about 50 meV, and thus, it will be possibly probed in planned neutrinoless double beta decay experiments. In addition, we mention the phases of the leptonic mixing matrix U . We find that the preferred value for the Dirac CP-violating phase δ is close to 2π in NO and close to π in IO, thus implying no leptonic CP violation in our scenario of the Zee model. The hint of δ ∼ 3π/2 from global fits to neutrino oscillation data, if confirmed, can therefore not be accommodated in any of the two orderings. We also find that the values of both Majorana CP-violating phases φ 1 and φ 2 are around 2π in both orderings.
As expected, the strong limits from other LFV processes imply that the NSI parameters χ m τ τ and ε m τ τ defined in eq. (3.36) are very suppressed in both orderings, i.e. χ m τ τ , ε m τ τ 10 −8 , and therefore well below future experimental sensitivity. The NSI parameters ε ρσ αβ , given in eq. (3.35), are also generated in the Zee-Babu model, where they reach values of about 10 −4 [89,102]. However, in the Zee model, NSIs turn out to be smaller, since neutrino masses are generated at one-loop level, while in the Zee-Babu model, the latter arise at two loops.
Furthermore, we check how our results change when imposing the fine-tuning parameter κ to be 10 instead of 1, see eq. (3.6). For both values of κ, the 1σ C.L. region is close to the upper limit on µ, even though the entire range down to µ = 1 GeV is allowed at 2σ C.L. The value of µ is not significantly correlated to the value of tan β. We find that the allowed ranges for the scalar masses (where the upper bounds determine the neutrino masses) depend critically on κ, having larger allowed mass ranges the larger the value of κ. Quantitatively, for κ = 1 at 1 (2)  Finally, we summarize some of our main results in table 7, where we display the 2σ C.L. regions for some of the most interesting observables and parameters for NO (for κ = 1, 10), and IO (for κ = 1, 10). We emphasize once more that for µ = 0 only upper bounds on the CLFV processes (and lower bounds on the scalar masses) exist, whereas for µ = 0, there are lower bounds for the CLFV processes (and upper bounds on the scalar masses). This means that the scalar sector cannot be arbitrarily heavy if neutrino masses are to be reproduced. The precise upper bound depends crucially on µ, see eq. (3.6). For µ = 0, the masses could be arbitrarily large, unlike the case of having µ = 0 and reproducing neutrino masses, which imposes that they are below about 2 TeV. We do not display the ranges for all observables, such as neutrino masses and leptonic mixing parameters, as their contributions to χ 2 min are shown in figure 2. In the Zee model (including µ = 0), the value of the muon AMM, which has not been included in the fit, is several orders of magnitude smaller than the experimental one. This implies a deviation of about 3.5σ. The allowed ranges for the scalar masses m A and m H are the same for both values of κ, but the two are not completely correlated (especially in IO). It is possible to simultaneously have m A 100 GeV and m H varying in the range 100 GeV m H 500 GeV. We conclude by stating that NO will be tested in the next generation CLFV searches, specially with τ → µγ and µe conversion, as well as searches for h → τ µ and the other new scalars at colliders.

Summary and conclusions
It is well known that there is LFV in the neutrino sector and this is also expected in the charged-lepton sector. In this work, we have studied the Zee model in detail, which is a simple extension of the SM that can accommodate neutrino masses and leptonic mixing if at the same time sizable signals in LFV processes are generated. We have performed a full numerical scan of the parameter space for three different cases (i) µ = 0, which implies massless neutrinos, (ii) NO, and (iii) IO. We have found that neutrino masses and leptonic mixing can be easily accommodated in NO, whereas IO is disfavored in comparison to NO due to the difficulty to fit the leptonic mixing angles θ 12 and θ 23 as well as the Dirac  CP-violating phase δ. In fact, if θ 23 turns out to be in the first octant, only NO would be allowed in the Zee model. Note also that none of the orderings can reproduce the hint of δ 3π/2 from global fits to neutrino oscillation data [60,61]. If expected sensitivities are achieved in τ → µγ and no signal is observed, a significant portion of the allowed parameter space for NO would be ruled out. This would put the Zee model under severe pressure, requiring an extension, e.g. involving the Yukawa couplings that give rise to terms proportional to m e and possibly large hierarchies among them, and relaxing the naturality demands on the trilinear coupling µ considerably, as Br(τ → µγ) ∼ 10 −9 would still rule out NO for κ ∼ 10. If no signals are observed in future µe conversion experiments, which are expected to increase their sensitivities by several orders of magnitude, the allowed regions in the parameter space will be strongly reduced, and IO will be basically excluded.
We have analyzed if the predicted rates of the Zee model for HLFV decays are observable at the LHC and future colliders. We have found that Br(h → τ µ) can be at the percent level, whereas Br(h → τ e) is at least two orders of magnitude smaller. If no signals are observed for Br(h → τ µ) at future colliders at the 10 −5 level, both orderings will be excluded at 1σ C.L. In general, we find that the expected sensitivities for CLFV will be more constraining in the near future. However, both CLFV and HLFV will have significant impact on the allowed parameter space of the Zee model.
In the model, neutrinoless double beta is due to only light neutrinos. Therefore, as usual, current experiments will be only sensitive to IO, while NO will only be tested if a further-order-of-magnitude improvement is achieved. In the Zee model, NSIs are always very suppressed, the strong limits from other CLFV processes and the fact that neutrino masses need to be generated at one-loop level.
In general, we have found that the masses of the new scalars should be at most a few TeV, which implies that they, especially the charged scalars that are pair-produced via Drell-Yan processes, can be searched for at the LHC. In particular, the masses of the neutral scalars and the charged scalar h + 1 are below 2.5 TeV, and typically they are lower than that, for both NO and IO. The phenomenology of the scalar sector is very modeldependent, like in general for 2HDMs, although it is possible to have sizable decays of the heavy neutral scalars into τ µ, correlated with the light Higgs ones.
We conclude by emphasizing that the general Zee model studied in this work is fully testable in the near future by combining different LFV processes. In particular, both orderings should be completely tested by CLFV and HLFV processes in the forthcoming years, as well as collider searches of the new scalars. Furthermore, if a signal in h → τ µ is observed at the LHC or in a future collider, the Zee model will be one of the best-motivated scenarios to accommodate it and at the same time describe neutrino masses and leptonic mixing.
In the Zee model, the parameter T is given by 13 where α em ≡ e 2 /(4π) is Sommerfeld's fine-structure constant 14 and the symmetric auxiliary function F is defined as Similarly, for the parameter S, following ref. [128] and adding the singly-charged contributions to the different gauge boson self-energies, we obtain and for the combination S + U , we find 15 where the renormalized auxiliary functions B 22 and B 0 are defined as ). 14 Note that παemv 2 = m 2 W s 2 W , where sW = sin θW and θW being the Weinberg angle. 15 We believe there are typos in the last two terms of ref. [128], which should have opposite signs.