Minimal 3-loop neutrino mass models and charged lepton flavor violation

We study charged lepton flavor violation for the three most popular 3-loop Majorana neutrino mass models. We call these models"minimal"since their particle content correspond to the minimal sets for which genuine 3-loop models can be constructed. In all the three minimal models the neutrino mass matrix is proportional to some powers of Standard Model lepton masses, providing additional suppression factors on top of the expected loop suppression. To correctly explain neutrino masses, therefore large Yukawa couplings are needed in these models. We calculate charged lepton flavor violating observables and find that the three minimal models survive the current constraints only in very narrow regions of their parameter spaces.


Introduction
One could understand the smallness of the observed active neutrino masses, in principle, if they are generated radiatively. It is therefore not surprising that loop models of neutrino masses have a rather long history [1][2][3][4]. Systematic classifications of loop models have been published for 1-loop [5], 2-loop [6,7] and, recently, even 3-loop [8] diagrams. For a detailed discussion we refer to the review [9].
In this work, we will study how upper limits on charged lepton flavor violating (CLFV) observables constrain 3-loop neutrino mass models. We will focus on some particular, wellknown models, which we consider "minimal" models. The term "minimal" here refers to the fact that for models at 3-loop level at least three different types of particles beyond the Standard Model (SM) particle content are needed, in order to avoid lower order diagrams. 1 The three models that we will study in this paper are the so-called cocktail [10], Krauss-Nasri-Trodden (KNT) [11] and Aoki-Kanemura-Seto (AKS) [12] models.
These three models are probably the best-known 3-loop models in the literature, and a number of other papers have studied them (or some variations thereof). The cocktail model, for example, has been studied also in [13]. There are also versions of the cocktail model in which the W bosons are replaced by scalars [14][15][16]. For the AKS model, one can find some discussion on phenomenology and vacuum stability constraints in [17][18][19], while a variant of the AKS model with doubly-charged vector-like fermions and a scalar doublet with hypercharge Y = 3/2 (plus the singlets of the AKS model) can be found in [20]. Other variants of the AKS model in which the exotic particles are all electroweak singlets can be found in [21,22]. Finally, for the KNT model, different phenomenological and theoretical aspects were studied in [23][24][25][26][27][28][29][30]. There are also variations of the KNT model, like the colored KNT [31][32][33][34], or a model with vector-like fermions added to the KNT model [35]. Other variants can be found in [36,37].
Common to all the three minimal models is that their neutrino mass diagrams are proportional to two powers of SM lepton masses. Together with the 3-loop suppression of 1/(16π 2 ) 3 , this results in the prediction of rather small neutrino mass eigenvalues, unless the new Yukawa couplings of the models take very large values. However, in all models off-diagonal entries for these new Yukawa couplings are required, since neutrino oscillation experiments have measured large neutrino angles, see for example [38] for a recent global fit of neutrino data. Therefore, one expects that CLFV limits will put severe constraints on these minimal models. This simple observation forms the motivation of the current paper.
The rest of this paper is organized as follows. In Section 2 we will set up the notation and briefly discuss two scalar extensions of the SM. In Section 3 we will discuss the cocktail model. We will first introduce the model and the neutrino mass generation mechanism in 3.1, and then we will present our numerical results for this model in 3.2. We start with the cocktail model, since the flavor structure of the neutrino mass matrix in this case is the simplest of the three models. We then discuss in a similar way the KNT model in Section 4 and the AKS model in Section 5. We close with a short discussion. A number of technical aspects on the calculation of the loop integrals are relegated to Appendix A.

Notation and conventions
In order to make the discussion more transparent for the reader, it is convenient to adopt a common notation and use the same conventions for the three models considered here. This is the aim of this section.
The three minimal 3-loop neutrino mass models studied in this work are based on the SM gauge group, SU(3) c × SU(2) L × U(1) Y . This local symmetry is supplemented by a global Z 2 parity, which is introduced to forbid the tree-, 1-and 2-loop contributions to the neutrino mass matrix, as explained below. The SM fermions as well as their charges under the SM gauge group are given in Table 1. They will all be assumed to be even under the global Z 2 symmetry. The particle spectrum of the 3-loop models explored in this paper may contain new fermions, and these will be fully specified for each model in the next sections. In what concerns their scalar sectors, they can be regarded as extensions of three well-known scenarios: the SM scalar sector, the Two Higgs Doublet Model (2HDM) scalar sector and the Inert Doublet Model (IDM) scalar sector. We now describe the scalar fields in these three minimal scenarios, the way the electroweak symmetry gets broken in each case and the Yukawa interactions with the SM fermions.
• SM scalar sector The SM scalar sector contains only the usual Higgs doublet, H, as shown in Table 2. This doublet can be decomposed in terms of its SU(2) L components as Table 3: 2HDM scalar sector, composed of the two SU(2) L scalar doublets Φ 1 and Φ 2 .
The Yukawa couplings of the SM are where we have defined H = iτ 2 H * , with τ 2 the second Pauli matrix. We have omitted flavor and SU(2) L indices in the previous expression to simplify the notation. The electroweak symmetry gets spontaneously broken by the Higgs vacuum expectation value (VEV), with v 246 GeV. In the three models discussed below, H will be even under the Z 2 parity.
• 2HDM scalar sector The 2HDM scalar sector is composed of two scalar doublets, Φ 1 and Φ 2 , with identical quantum numbers under the SM gauge symmetry, as shown in Table 3. They can be decomposed in terms of their SU(2) L components as Since both scalar doublets have exactly the same quantum numbers, and in particular since they will both be assumed to be even under the Z 2 symmetry, flavor changing neutral current interactions are in principle present. This dangerous feature can be fixed by introducing a second (softly broken) Z 2 symmetry, under which one of the two doublets and some of the SM fermions are charged. There are several possibilities, and here we will just assume that this symmetry makes Φ 1 leptophilic, and Φ 2 leptophobic. 2 Under this assumption, the 2HDM Yukawa interactions are given by Again, flavor and SU(2) L indices have been omitted for the sake of clarity. We see that, as explained above, Φ 1 only couples to leptons, while Φ 2 only couples to quarks. In the 2HDM, both scalar doublets are assumed to take VEVs, such that the usual electroweak VEV v is given by We also define the ratio • IDM scalar sector In the IDM, a second scalar doublet denoted as η is introduced. In contrast to the 2HDM, this doublet is odd under the Z 2 parity, as shown in Table 4. The inert doublet η can be decomposed in terms of its SU(2) L components as Since the SM fermions are even under Z 2 , η does not couple to them and the IDM Yukawa interactions are exactly the same as those in the SM, see Eq.
(2). The scalar potential of the IDM is assumed to be such that only the SM Higgs doublet takes a VEV, Therefore, electroweak symmetry breaking takes place in the standard way and the Z 2 parity remains exactly conserved. Finally, one can split the neutral component of the η doublet as so that η R and η I are, respectively, the real and imaginary parts of η 0 . Under the assumption of CP conservation in the scalar sector, these two states are mass eigenstates, since the Z 2 symmetry forbids their mixing with the SM neutral scalar.
In the next sections we will completely specify the new fields and interactions for the three models considered here. In what concerns the leptonic sector, some comments are in order. Neutrino oscillation data fixes the mass squared splittings ∆m 2 Atm and ∆m 2 , the three leptonic mixing angles and the so-called "Dirac" phase δ. We will work in the basis in which the charged lepton mass matrix is diagonal, This implies that the unitary matrix U that brings the 3 × 3 Majorana neutrino mass matrix M ν to diagonal form as corresponds to the leptonic mixing matrix that enters the charged current interactions. We will use the PDG parametrization [39] and write U as where R ij are the standard rotation matrices and P is a diagonal matrix containing the Majorana phases, P = diag(1, e iα 12 /2 , e iα 13 /2 ) .

Cocktail model
We begin with the so-called cocktail model, introduced in [10], since the neutrino mass matrix in this model has the simplest flavor structure.

The model
The cocktail model can be regarded as an extension of the IDM. In addition to the IDM fields, the particle content of the cocktail model includes the two SU(2) L singlet scalars S and ρ, singly and doubly charged. Interestingly, the model does not have any new fermion, just scalars. The η and S scalar fields are taken to be odd under the Z 2 parity, while the rest of the fields in the model are even. 3 The quantum numbers S and ρ are given in Table 5.
Counting also η, there are then three new multiplets in the cocktail model, with respect to the SM. The Lagrangian of the cocktail model contains only one additional Yukawa term with respect to the SM, where h is a symmetric 3 × 3 matrix. Flavor indices have been omitted in this expression for the sake of clarity. In addition, the new scalar potential couplings are given by We have omitted SU(2) L indices to simplify the notation. The parameters µ 1 and µ 2 are trilinear couplings with dimensions of mass, κ and all λ's are dimensionless. Most important is the term proportional to λ 5 , see discussion below. The singly charged scalars S + and η + mix after electroweak symmetry breaking, due to the term proportional to µ 1 . This leads to two H + i mass eigenstates, with mixing angle β, see Appendix A.1. The model also includes the doubly-charged scalar ρ ++ , with mass The cocktail model has other interesting features that will not be discussed in any detail here. For instance, the Z 2 parity of the model is conserved after electroweak symmetry breaking, so that the lightest Z 2 -odd state is stable and can in principle constitute a good DM candidate.

Neutrino masses
The 3-loop diagram leading to neutrino masses in the cocktail model is shown in Fig. 1. In the unitary gauge this diagram is the only diagram contributing to the neutrino mass matrix. However, in order to understand how to maximize the contribution of this diagram to the neutrino mass matrix, it is more useful to calculate all diagrams in Feynman-'t Hooft gauge. This is discussed in detail in Appendix A.1.
In an analogous way to the well-known scotogenic model [40], the diagram shown in Fig. 1 vanishes in the limit m 2 η R − m 2 η I ∝ λ 5 → 0, since in this limit the model conserves lepton number. We can then write the neutrino mass matrix in the cocktail model as:  where m i and m j are charged lepton masses. Here we have hidden all the complexities of the calculation in the dimensionless factor F Cocktail . This factor contains the loop integrals, depending on the masses of the scalars, and prefactors containing coupling constants, etc, see Appendix A.1.

Results
The cocktail model is an example of a type-II-seesaw-like model. In this class of models, the neutrino mass matrix is proportional to a symmetric Yukawa matrix, with Y ij = Y ji and Λ some generic mass scale. This allows one to fit the observed neutrino masses and mixing angles in a trivial way. Furthermore, the tight relation given in Eq. (20) implies very specific predictions for ratios of CLFV observables and strongly reduces the number of free parameters in the model. As we will discuss now, this has very important consequences.
From the experimental data we can reconstruct the neutrino mass matrix in the flavor basis as This allows us to calculate the Yukawa h necessary to fit the experimental data using the expression in Eq. (19). We find  Figure 2: Yukawa couplings h ij as function of the lightest neutrino mass, calculated with F max Cocktail . These Yukawas should therefore be understood as lower limits. In both plots we have used the best fit point data from the global oscillation fit [38], except δ = 0. The plot to the left shows the case (α 12 , α 13 ) = (0, 0), the plot on the right (α 12 , α 13 ) = (π, 0). The dashed gray (black) lines in the background are rough estimates for the typical size that the Yukawa couplings should have, in order to satisfy limits from muon (tau) CLFV decays. These lines are only for orientation.
where M e is the diagonal matrix with the measured charged lepton masses, see Eq. (12).
Let us first make a rough numerical estimate. Choosing normal hierarchy (m ν 1 → 0) and δ = 0 for simplicity and inserting m ρ ++ = 800 GeV, which is roughly the current experimental bound from LHC data [41][42][43] These values are obviously much too large to be realistic. We therefore searched the parameter space, intending to identify regions, in which h can fulfill the bounds from perturbativity and lepton flavor violation searches. This search was done in two steps. First, we maximize F Cocktail and λ 5 . For λ 5 we use λ 5 = 4π, the largest value allowed by perturbativity. We then scanned all free mass parameters entering in F Cocktail , for details see appendix A.1. Generally speaking, F Cocktail is maximized when µ 1 , µ 2 and κ take the largest values allowed, while the remaining free mass eigenvalues of the model take the lowest possible values allowed by experimental searches. The maximal value of F Cocktail found in this numerical scan is F max Cocktail 192. We will use this number in all plots below. This choice is conservative in the sense that the Yukawa couplings h ij will be larger for all other choices, thus constraints from charged lepton flavor violation searches will only be more stringent in other parts of the parameter space.
Once F max Cocktail is fixed, one can scan over the free parameters in the neutrino sector. Oscillation data [38] fixes rather well ∆m 2 Atm , ∆m 2 and all three mixing angles; there is also an indication for a non-zero value of δ. This leaves us with three essentially free parameters, the two Majorana phases and m ν 1 , equivalent to the overall neutrino mass scale, for which there are only upper limits from neutrinoless double beta (0νββ) decay [44,45] and cosmology [46]. We note that there is a slight preference in the data for normal hierarchy (NH, also called normal ordering) over inverted hierarchy (IH).
In Fig. 2 we plot the absolute values of the 6 independent entries in the Yukawa matrix h as a function of m ν 1 . The oscillation data have been fixed at their best fit point (b.f.p.) values, except δ = 0 for simplicity. The plot on the left was obtained with vanishing Majorana phases, whereas the one on the right takes (α 12 , α 13 ) = (π, 0). Given that we used F max Cocktail in this plot, the numerical values of h ij are much smaller than in Eq. (23), but h 11 is still in the non-perturbative region everywhere in the left plot. In the right plot, however, there are two special points, where cancellations among different contributions of the neutrino mass eigenstates lead to a vanishing value for either h 11 or h 12 . Such cancellations are well known in studies of 0νββ decay. The effective Majorana mass, m ee , 4 depends on the Majorana phases in the same way as h 11 . As in m ee , one can therefore not obtain a cancellation for the cases (i) NH without Majorana phases, and (ii) IH for any choice of parameters. The cocktail model can therefore explain neutrino data only for normal hierarchy and some particular combination of Majorana phases, as we are going to discuss now in some more detail.
As Fig. 2 demonstrates, only in some exceptional points can h 11 be small enough to enter the perturbative region. We therefore scanned over (α 12 , α 13 ) and m ν 1 , in the full 3 σ range of oscillation data. In this scan, we calculate the CLFV observable Br(l i → l j l k l m ), with different combinations of lepton flavors, for the minimal value of m ρ ++ allowed by LHC data. Fig. 3 to the left shows Br(l i → l j l k l m ) for all different combinations of i, j, k, m using the b.f.p. of neutrino oscillation data. The most stringent constraint on the model comes from the experimental upper limit on Br(µ → 3 e)≤ 10 −12 [47]. The plot to the right then shows the allowed regions in parameter space, scanning over the complete range of oscillation parameters and phases. All acceptable points lie in the range m ν 1 = (2 − 10) meV. Fig. 4 shows a scan over the allowed range of Majorana phases and the lightest neutrino mass. The plot to the left shows the plane (α 13 , α 12 ), the one to the right (α 12 , m ν 1 ). The model can fulfill the constraint from Br(µ → 3 e) only in a very narrow range of phases. In Figure 4: Allowed parameter space for α 12 , α 13 and m ν 1 . Note that m ν 1 is shown in units of meV. Neutrino oscillation data was scanned over the 3 σ uncertainties, except δ which is taken at its best fit value for simplicity. particular, α 12 has to be close to π in all points, while also m ν 1 is fixed in a rather narrow interval.
Finally we note that the acceptable points of the model lie in regions of parameter space where the 0νββ decay observable m ee is unmeasurably small. There is, however, a 1-loop short-range diagram contributing to 0νββ decay in the cocktail model [48], see Fig. 5. This diagram depends on the same parameters as the 3-loop neutrino mass diagram in Fig. 1.
In particular, note that the sum over η R,I generates the same dependence on λ 5 as for the neutrino mass.
We have calculated this diagram and estimated its contribution to the 0νββ decay halflife, including the QCD running of the short-range operator [49,50]. Using the same mass parameters that maximize the 3-loop diagram, in particular m ρ ++ = 800 GeV, the current limit on the half-life of 136 Xe [44] imposes a limit on h 11 of roughly |h 11 | < ∼ 5 × 10 −4 . This limit is around a factor ∼ 7 more stringent than the one obtained from the upper limit on Br(µ → 3 e). 5 We can conclude that the cocktail model is severely constrained from perturbativity arguments and from searches for CLVF. The model has acceptable points only within a narrow window of m ν 1 and for particular combinations of the Majorana phases.

KNT model
We continue with the KNT model [11]. This was the first radiative neutrino mass model at 3-loop order proposed.

The model
In addition to the SM particles, the KNT model contains three copies of the fermionic singlet N and two singly-charged singlet scalars X and S. A discrete Z 2 symmetry is imposed, under  which S and N are odd and the rest of particles in the model are even. The quantum numbers of the new particles in the KNT model are given in Table 6.
The Lagrangian of the model contains the following pieces where we have omitted SU(2) L and flavor indices to simplify the notation. We note that f is an antisymmetric 3 × 3 Yukawa matrix, while M N is a symmetric 3 × 3 Majorana mass matrix, which we take to be diagonal without loss of generality. The scalar potential of the model also contains additional terms besides those in the SM. These are given by The presence of the λ S quartic coupling precludes the definition of a conserved lepton number. Indeed, one can easily see that the simultaneous presence of the Lagrangian terms in Eqs. (24) and (25) breaks lepton number in two units. The masses of the physical scalar states in the KNT model are given by We also note that the lightest Z 2 -odd state in the KNT model is completely stable. Assuming the hierarchy M N 1 < m s 1 < m s 2 , this state is the lightest fermion singlet, which then constitutes a good DM candidate. In fact, the KNT model is historically the first radiative neutrino mass theory with a stable DM candidate running in the loop.

Neutrino masses
The Z 2 symmetry forbids the standard Higgs Yukawa coupling with the lepton doublet L and the N singlets. Therefore, the usual type-I seesaw contribution at tree-level is absent. Instead, neutrino masses are generated at 3-loop order as shown in Fig. 6. The neutrino mass matrix is given by where m α is the mass of the α charged lepton and F KNT is a loop function that depends on the masses of the scalars and fermions running in the loops. More information about this function can be found in Appendix A.2.
It is important to stress that in the KNT model, each entry in (M ν ) ij contains the sum over the SM charged lepton masses. Therefore, different from the other models discussed in this paper, the suppression of the entries in M ν is at most m 2 µ . The neutrino fit can then reproduce experimental data with Yukawas which are considerably smaller than in the cocktail or AKS models.

Results
We start this section again with a discussion of the neutrino mass fit. The coupling f in Eq. (24) is antisymmetric, thus the determinant of the neutrino mass matrix in Eq. (29) is zero, implying that one neutrino is massless. This is reminiscent of the 2-loop Babu-Zee model of neutrino mass [3,4], where the same singly charged scalar is used. In our fitting procedure we use therefore an adapted version of the solution found in [51,52] for the Babu-Zee model.
The procedure consists of two steps. First, because det(f ) = 0, the matrix has one eigenvector a = (f 23 , −f 13 , f 12 ), which is also an eigenvector of M ν : This implies three equations, one of which is trivial, while the other two allow to express the ratios (f 13 /f 12 , f 23 /f 12 ) as functions of the neutrino angles and phases only. These solutions depend on the neutrino mass hierarchy. Next, we can write the neutrino mass matrix as c contains all global constants, we have used f T = −f and M aux is an auxiliary matrix, which is complex symmetric. This defines a set of 6 complex equations relating the entries in M aux to neutrino data. With three independent entries f ij , we can use three of the six equations to express three entries in M aux as a function of the remaining ones, neutrino data and f ij . The resulting equations are very lengthy and not at all illuminating, so we do not present them here. The definition of M aux in Eq. (31) shows that where and With M aux being complex symmetric, we can use a suitably modified [53,54] Casas-Ibarra parametrization [55] to express the matrix g as M aux and U aux are the eigenvalues and eigenvectors of the auxiliary matrix M aux . Although in principle it would be possible to determineM aux and U aux in terms of the input neutrino data analytically, in practice we find these two matrices numerically for any input point of experimental data and choice of free parameters. Finally, R is a 3 × 3 orthogonal matrix.
In summary, neutrino oscillation data provides 6 constraints: ∆m 2 Atm , ∆m 2 , three angles and the CP-phase δ. A number of free parameters can then be scanned over, using the above procedure. In the neutrino sector we still have α 12 . 6 The matrix f is fixed from experimental data, up to the overall scale of the matrix. We choose f 12 as the free parameter. The matrix R, in the most general case, contains 3 complex angles. There are 3 right-handed neutrino masses, M N i , and 2 scalar masses, m s 1,2 . And, finally, we can use 3 of the 6 equations for M aux to eliminate some particularly chosen (M aux ) ij . This leaves as free inputs the remaining 3 entries in M aux .
Up to now, we have been completely general in our discussion. However, there is still a certain freedom as to which 3 entries in (M aux ) ij we fix via 3 of the equations defined by Eq. (31). In practice, we choose to solve for (M aux ) 22 , (M aux ) 23 and (M aux ) 33 and assume (M aux ) 1k = 0. This particular choice is motivated by the observation that in this limit all terms in g ∝ 1/m e disappear. In other words, this solution guarantees that the contribution to µ → eγ and τ → eγ from loops involving s 2 and N i are automatically absent in our scans, due to g 1k = 0 ∀k. Our ansatz is therefore the optimal choice for minimizing fine-tuning on the other parameters in g. 7 Before we explore the remaining parameter space of the model, we must consider lower limits on the masses of the charged scalars from accelerator searches. LEP provides a lower limit on charged particles decaying to leptons plus missing momentum, which will essentially rule out all values of m s 1 below 100 GeV [39] and similarly for m s 2 , unless M N 1 is close to m s 2 . 8 At the LHC there is currently no specific search for particles with the quantum numbers of S and X. However, slepton pair production with their subsequent decays to a lepton plus a neutralino provide the same signal and thus, we can make a reinterpretation of the corresponding searches at CMS [56] and ATLAS [57]. The CMS slepton search [56] is based on 35.9/fb, while ATLAS' chargino and slepton search [57] uses 139/fb. The ATLAS limits are correspondingly more stringent and we will therefore discuss these. We have implemented the KNT model in SARAH [58,59] and generated SPheno routines [60,61] and model files for MadGraph [62][63][64]. We have then calculated cross sections with MadGraph to recast the results of [57]. For s 1 the mass range between (very roughly) m s 1 = (250 − 400) GeV is excluded by this search. The range m s 1 = (100 − 250) is currently unconstrained, due to large backgrounds in [57]. For s 2 the limits are even weaker, unless m s 2 − M N 1 is larger than (50 − 70) GeV, depending on m s 2 .
Let us now turn to the discussion of CLFV. Consider first the antisymmetric Yukawa coupling f . Neutrino data requires all three elements of f to be non-zero, thus there will always be a non-zero value for the three possible decays µ → eγ, τ → eγ and τ → µγ. The constraint from µ → eγ is the most stringent one. However, since we still have the overall scale of f as a free parameter, in our choice the value of f 12 , we can use it to fix Br(µ → eγ) to 6 Since one neutrino is massless, only one of the two Majorana phases, i.e. α 12 , is physical. 7 We have explored other ansätze but concluded that this is indeed the optimal one. For instance, we can generate textures with either the 2nd or 3rd column of g vanishing. These will make the τ → eγ or τ → µγ branching ratio vanish, but at the cost of a large µ → eγ branching ratio. We have also considered a solution with any or all of (M aux ) 1k = 0. However, in this case we did not find any configuration for the remaining free parameters that induces a cancellation that suppresses the µ → eγ branching ratio. 8 There are no accelerator limits on N , since the Z 2 symmetry prohibits their mixing with the active neutrinos. To the left, with current experimental limits; to the right for the expected future experimental limits. The required size of the Yukawas is color-coded. Bluish-grey points mean that at least one entry in g is larger than 4π.
the upper limit (present or future) for any point in the parameter space. Since the neutrino mass matrix is proportional to the square of the matrices f and g however, once this choice is made there is no longer any overall scaling freedom in the coupling g. Putting the calculated Br(µ → eγ) to equal the experimental bound will generate the smallest values for the entries of g allowed in the model parameter space. A smaller upper limit on Br(µ → eγ) will lead to larger g and thus more stringent constraints from τ → µγ.
We then scanned over the remaining parameters of the model numerically. Consider first the case of NH. Some examples for Br(τ → µγ) are shown in Fig. 7. In this plot we have chosen the fixed value m s 1,2 = 100 GeV, corresponding to the experimental lower limit, and the three right-handed neutrino masses all equal to a common M N . 9 The points are scanned over the allowed 3 σ ranges for the neutrino data for NH. The size of the largest entry in g is color-coded in the points. The plot to the left has been calculated for the current experimental limit on Br(µ → eγ)< 4.2 × 10 −13 [65], the plot to the right is for the expected future limit Br(µ → eγ)< 6 × 10 −14 [66]. For the choice of m s 1,2 = 100 GeV no valid point with g ij ≤ 4π ∀ij remains in the parameter space. Constraints are more stringent for IH, and therefore the same conclusion is reached.
We therefore scanned over m s 1,2 ≡ m S and M N i simultaneously. The results are shown in Fig. 8. Here, M N i are varied within 20% of a common M N . The range of M N is color-coded in the points. Again, the plot to the left is for the current bound on Br(µ → eγ), while the plot to the right is for the future bound. In these plots, points with non-perturbative couplings are shown in bluish color. This bound eliminates all points below roughly M N = O(100) GeV already with the current experimental bound on Br(µ → eγ), see however the discussion below. We show only the cases with a trivial R matrix. For non-zero angles in R the results look similar, although fewer points lie in the perturbative regime.
The combined constraints of perturbativity and future limits from CLFV searches would put a lower bound on m S roughly of order (180 − 200) GeV. This limit becomes stronger for lower values of M N , as the plots shows.   The above discussion is strictly valid only for the case where the three right-handed neutrinos have similar masses. For hierarchical right-handed neutrinos the constraints are usually dominated by the lightest of these. There exist, however, exceptional points in the parameter space, where the contributions to Br(µ → eγ) from the three different neutrinos conspire to (nearly) cancel each other. This is shown in Fig. 9. The figure shows Br(τ → µγ) as a function of the "lightest" right-handed neutrino mass, for different choices of M N 2,3 . Br(τ → µγ) is dominated by the lightest mass eigenstate, except in some particular points, where cancellations occur. Figs. 7 and 8 do not cover these exceptional combinations of parameters.
We have repeated the scans discussed above also for the case of IH. An example is shown in Fig. 10. IH requires larger Yukawa couplings, since now two neutrino have masses of order ∆m 2 Atm . Thus, many more points in the parameter space are ruled out due to the perturbativity constraint. This pushes both, fermion as well as the scalar, masses to larger values. Indeed, already with current constraints there are no points with m S below roughly 600 GeV.
Finally, let us mention that in the KNT model there is no short-range diagram contribut- Figure 10: Same as Fig. 8 but for IH for the neutrino masses. Bluish points are excluded due to perturbativity arguments.

AKS model
A general class of models is represented by the AKS model [12]. In this case the particle content is extended to include new scalars and fermions.

The model
The AKS model extends the usual 2HDM with the real scalar singlet ϕ, the singly charged scalar S and three generations of singlet fermions N . Even though a more minimal version with only two generations of N is possible, we will consider three in the following. The fields S, ϕ and N are assumed to be odd under the Z 2 parity, while the rest of the particles are even. The quantum numbers of the new particles in the AKS model are given in Table 7.
As explained in Sec. 2, an additional softly-broken Z 2 symmetry is introduced to avoid dangerous flavor changing neutral currents. We choose to follow [12] and use this symmetry to couple one of the scalar doublets (Φ 1 ) only to leptons, and the other (Φ 2 ) only to quarks.
Due to this choice, the Yukawa couplings of the model are given in Eq. (5), along with the Yukawa −L ⊃ Y * N c e R S + h.c. .
One can also write Majorana masses for the N singlets, with M N a symmetric matrix. The scalar potential of the model is given by As usual, we have omitted SU(2) L indices in the previous expression. We point out that lepton number would be restored in the limit κ → 0. The presence of this coupling breaks lepton number in one unit. After electroweak symmetry breaking, the doublet scalars Φ 1 and Φ 2 get mixed. The mass eigenstates resulting from this mixing are the SM Higgs, another Higgs, a new charged scalar, and a pseudoscalar. The charged and neutral Goldstone bosons are absorbed by the Z and W gauge bosons. The mass matrix for the CP-even neutral states in the basis H 0 = Re (Φ 0 1 , Φ 0 2 ) T is given by The CP-odd neutral scalar mass matrix in the basis One finds a massless state, the Goldstone boson that becomes the longitudinal component of the Z boson. The other state has a mass while the mass of the Z boson is m 2 Z = 1 4 v 2 (g 2 1 + g 2 2 ). The mass matrix for the charged states in the H ± = (Φ ± 1 , Φ ± 2 ) T basis is  Again, after diagonalization one obtains a massless state, identified with the Goldstone boson that becomes the longitudinal part of the W boson, and a massive physical charged scalar with mass The mass of the W boson is given by the standard expression m 2 W ± = 1 4 g 2 2 v 2 . Finally, the masses of the singlet scalars ϕ and S are As in the cocktail model, the lightest Z 2 -odd state in the AKS model is stable and can constitute a DM candidate.

Neutrino masses
In the AKS model, neutrino masses are induced at 3-loop order, as shown in the diagrams of Fig. 11. The resulting neutrino mass matrix is given by where m i is the mass of the i-th charged lepton and F AKS is a dimensionless loop function that depends on the masses of the scalars and fermions in the loop. More details about the calculation of this loop function can be found in Appendix A.3. The Yukawa matrix Y in the AKS model does not have any specific symmetry. Therefore, this model represents the general class of models in which the Yukawa matrices can be described by using a generalization of the Casas-Ibarra parametrization [55] (see also [53,54]).

Results
The Yukawa structure of the neutrino mass matrix shown in Eq. (45) resembles that of the type-I seesaw. In order to fit the experimental oscillation data, we use the Casas-Ibarra parametrization introducing the neutrino mass matrix in the flavor basis given in Eq. (21). We find that where M N has been taken to be diagonal and R is an arbitrary complex 3 × 3 orthogonal matrix. We include F AKS as it is, in general, a function of the eigenvalues of M N , see clearly in the non-perturbative regime. Insisting on perturbative Yukawa couplings thus calls for cancellations, especially in the first column, proportional to 1/m e . Moreover, even if for a choice of parameters, the Yukawa lives at the edge of perturbativity, one should take care of the constraints coming from CLFV. Especially µ → eγ, given the hierarchy among the entries of the Yukawa matrix Y . In order to avoid non-perturbativity and CLFV constraints, first we exploit the freedom in R. We fix two of the complex angles to make two entries of the Yukawa matrix zero or close to zero. We choose Y 21 and Y 31 . With this, we find that the third free angle of R is not enough to cancel another entry in the Yukawa matrix. Therefore, we can only fix the values of the phases and m ν 1 to minimize or cancel Y 11 , similarly to the cocktail model, or Y 12 , to live below the experimental limit on Br(µ → eγ), proportional to |Y k1 Y * k2 | 2 . From now on, we also consider κ = 4π, at the edge of perturbativity, and tan β = 1.
In Fig. 12 to the left, we show the behavior of the first row of Y for α 12 = α 13 = δ = π. We considered for simplicity that all the scalar masses are equal to m S + = m ϕ = m H ± ≡ m S and all the N singlet fermion masses to be degenerate, m N , and minimize m N /F AKS to find the lowest value of Y , see Eq. (46). We found this minimum for m N = 272 GeV and m S = 100 GeV, where F AKS ≈ 0.44, compatible with the limit on scalar masses from LEP [39]. Here we do not show the other four non-zero Yukawas for simplicity. They are nearly constant and of order 0.1. Similar to the cocktail model (Fig. 2), poles exist in the different Yukawa entries for particular values of the phases and m ν 1 . The main difference lies in the divergence that appears when Y 11 = 0. This is caused by our choice of R matrix, such that Y 21 = Y 31 = 0. In this case, the pole in Y 11 does not imply a pole in Br(µ → eγ) or Br(τ → eγ), as it can be seen in Fig. 12 to the right. In fact, the product of |Y 11 Y * 13 | remains constant over the pole and very close to the current experimental limit of 3.3 × 10 −8 [67]. Only the region around the pole in Y 12 is allowed by the experimental limit Br(µ → eγ)< 4.2 × 10 −13 [65].
To sum up, the parameter space of the AKS model is constrained mainly by perturbativity and Br(µ → eγ). The former can be addressed with the freedom in R to set Y 21 and Y 31 to  Figure 12: To the left, the entries of the first row of the Yukawa coupling matrix Y zoomed around the poles for Y 11 and Y 12 for α 12 = α 13 = δ = π. To the right, the calculated Br(l α → l β γ). Both computed fixing R for Y 21 = Y 31 = 0 and at the minimum allowed value of m N /F AKS . While a pole for Y 11 exists, no pole for Br(µ → eγ) or Br(τ → eγ) is associated to it, due to the divergence of Y 12 and Y 13 on the pole. Note that only near the pole for Y 12 Br(µ → eγ) is below the experimental limit.
zero. As well as by fixing the Majorana and Dirac phases, and the lightest neutrino mass, to be near the pole of Y 11 , where its value is lower than 4π. On the other hand, to be below the experimental limit on Br(µ → eγ), a similar fine-tuning of the phases and m ν 1 should be done to be around the narrow pole of Y 12 . The parameter space is then restricted to those values of the phases and m ν 1 where the poles of Y 11 and Y 12 exist, and they are close enough to each other to avoid the limit on Br(µ → eγ) while Y 11 is still perturbative. In Fig. 13 we show the value of Br(µ → eγ) scanning over the complete range of oscillation parameters (NH) and phases. On the right, we give the limit due to perturbativity of Y 11 , reducing the parameter space to a small window of m ν 1 = (4.5 − 20) meV. Note that like in the cocktail model, Y 11 behaves as m ee , and for m ν 1 > ∼ 10 meV, m ee has no pole, so Y 11 is in the non-perturbative region. Moreover, the cancellation of Y 11 and Y 12 only occurs for NH, so the model can only explain neutrino data with this neutrino mass ordering. In the following, we shall consider only NH. Fig. 13 not only implies a constraint on m ν 1 , but also on the phases. In Fig. 14 we show the points allowed by perturbativity and the experimental limits on Br(µ → eγ), Br(τ → eγ) and Br(τ → µγ), for the values of the three phases. We scanned over the phases and masses, with m S > 100 GeV, allowing oscillation data to vary in 3σ. As it can be seen, α 12 should be closely around π, while δ is constrained to values between roughly π/2 and 3π/2. For δ outside this window, there is no cancellation of Y 12 . In the following, we restrict the results shown to the region where Br(µ → eγ)< 4.2 × 10 −13 . Now we move to analyse τ → eγ and τ → µγ. As shown in Fig. 12 (right), while Br(τ → µγ) is below the experimental limit, except on the pole of Y 11 , Br(τ → eγ) is mainly constant and close to the experimental limit. In Fig. 15 we give both branching ratios fixing δ to the b.f.p. and scanning over the uncertainties in the rest of the oscillation parameters. We consider (m N /F AKS ) min with m S = 100 GeV and m N = 272 GeV. Points colored in gray correspond to non-perturbative Yukawas. We see that while Br(τ → µγ) is safe, the allowed region on the left plot is severely constrained by the experimental limit Br(τ → eγ)< 3.3×10 −8 . This tension can be mitigated by raising the masses, see Fig. 16. For 10 -2 Figure 13: Br(µ → eγ) scanned over neutrino oscillation data in 1σ (light blue) and 3σ (dark blue) ranges. This plot scans over the Majorana phases. The shaded gray area corresponds to the most conservative limit to non-perturbative Yukawas. Figure 14: Allowed parameter space for α 12 , α 13 and δ. Neutrino oscillation data was scanned over the 3σ uncertainties, except δ which was left free. For points outside this region, either Y 11 is non-perturbative or Br(µ → eγ) is above the experimental limit.  Here, we are maximizing the allowed parameter space in terms of α 13 . For the gray points, at least one entry in Y is larger than 4π.
the AKS model, the dominant contribution to Br(l α → l β γ) is approximately proportional to 1/M 4 , with M the dominant scale [68]. On the other hand, m N /F AKS is minimal for masses around m S = 100 GeV and m N = 272 GeV. So for masses away from these values, m N /F AKS increases and, consequently, the absolute scale of the Yukawas increases as well (see Eq. (46)), hence narrowing the region where the Yukawas are perturbative. For m N (m S + ) ∼ 10 6 GeV, we found no points allowed by perturbativity and the experimental limit on Br(µ → eγ). In Fig. 16, in order to minimize the Yukawas, we fixed m ϕ = m H ± = 100 GeV and change m S + and m N , which enter in the calculation of Br(l α → l β γ).
A similar analysis can be done scanning over the Majorana phases too. Fig. 17 shows Br(τ → eγ) as a function of α 13 for different fermion and scalar masses. The allowed parameter space is bigger for m N around 272 GeV, where m N /F AKS is minimal. For different masses the parameter space narrows, because m N /F AKS increases, as explained before. The upper limit is due to the phenomenological limit m S + > 100 GeV, as for m N m S + , Br(τ → eγ) is dominated by m S + . On the other side, while going to larger m N reduces considerably Br(τ → eγ), a lower limit always exists due to perturbativity. We close this discussion with a short comment on 0νββ decay. There is no short-range diagram for 0νββ decay in the AKS model. Since, as discussed above, the AKS model survives only for normal hierarchy and in the part of parameter space where m ee is largely cancelled, observation of 0νββ decay in the next round of experiments would definitely rule out AKS as an explanation of neutrino masses.

Discussion
In this paper we have considered the cocktail, KNT and AKS models and studied their CLFV phenomenology. In these models, Majorana neutrino masses are generated at the 3-loop order, which naturally implies that large Yukawa couplings are required in order to reproduce the mass scales observed in neutrino oscillation experiments. As a result of this, perturbativity is typically lost. We have shown that one can decrease the Yukawa couplings by tuning some of the free parameters of these scenarios, such as the lightest neutrino mass m ν 1 or the Dirac and Majorana phases contained in the leptonic mixing matrix U . However, even after these parameters are tuned to recover perturbativity, the resulting CLFV branching ratios tend to largely exceed the existing bounds. In order to reduce the CLFV rates further tuning is needed. Our main conclusion is that the three models survive only in tiny correlated regions of their parameter spaces.
One should note that CLFV alone cannot exclude any of these models. The reason is that one can always reduce the CLFV rates as much as necessary by tuning the parameters of the model more finely. However, additional experimental handles exist. First, perturbativity imposes upper limits on the masses of some of the particles running in the loops. The reason is simple: larger mediator masses would imply a stronger suppression of the loop functions and then require larger Yukawa couplings. Thus, also future searches at the LHC in the high-luminosity phase would further restrict the available parameter space. An important experimental handle on the models is 0νββ decay. Since m ν 1 and the Majorana phases must be tuned for the models to survive, the effective 0νββ neutrino mass m ee becomes strongly constrained and definite predictions for the 0νββ rates are obtained for the AKS and KNT models. Any observation of 0νββ decay with the next generation of experiments would definitely rule out the AKS model. For the KNT model, because m ν 1 0, m ee has to be either in the range m ee (2 − 6) meV or (15 − 50) meV for normal hierarchy or inverted hierarchy. Only the cocktail model is more flexible in its predictions for 0νββ decay, due to additional contributions from a sizable short-range diagram.
We mention also that only the KNT model can explain neutrino data for both hierarchies. Neither the cocktail nor the AKS model has any acceptable point in all of their parameter space in the case of inverse hierarchy.
A crucial ingredient in our analysis is the allowed size for the quartic scalar potential couplings that play a role in the neutrino mass generation mechanism, for example λ 5 in the cocktail model, λ S in the KNT model and κ in the AKS model. Since neutrino masses are proportional to (some power of) these couplings, the larger they are, the smaller the Yukawa couplings can be. In our analysis, scalar couplings as large as 4π have been allowed. A more restrictive choice, with couplings at most of O(1), would alter the conclusions dramatically. In fact, all three models would already be ruled out, if all their couplings are restricted to be not larger than O(1).
Finally, we emphasize again that our strong claims only apply to the three minimal models considered here. There are several ways to modify these models so that they can evade the perturbativity and flavor constraints. For instance, one can introduce new exotic states in order to get rid of the proportionality to the charged lepton masses, at the origin of the problems discussed in our paper. Also, one may enhance the contributions to the neutrino mass matrix by using colored states. Nevertheless, we also note that there may be many other 3-loop (or 4-loop) neutrino mass models with the same issues.

A Loop integrals
In this Appendix we discuss the calculation of the loop integrals in the cocktail, KNT and AKS models. Here, we derive the loop functions used in the previous sections. For their computation, we did not rely on approximations, but implemented the full integral numerically using pySecDec [69].
Moreover, all the integrals shown here, can be factorized in terms of five master integrals (see for example [70]), as normally done. Nevertheless, we decided not to do it, because there is still no analytical general solution to all the 3-loop master integrals and their factorization could lead to numerical precision issues. Note that these five master integrals have divergent parts, while the full integral is finite.

A.1 Cocktail model
To compute the dimensionless integral F Cocktail in Eq. (19), we choose the Feynman-'t Hooft gauge ξ = 1. In this gauge, the propagator of the W µ boson has no momenta structure in the numerator, while the standard Goldstone H + contribution with a mass m 2 W should be  Figure 18: Dimension 5 mass diagrams in the gauge basis. Note that given the chirality, the corresponding integral has two momenta in the numerator. This is denoted with / ∂. When referring to these diagrams we will use the notation I included. We decided to show the diagrams in the gauge basis to be able to identify the different contributions that enter in F Cocktail .
We identified 12 different diagrams in the gauge basis with dimensions 5, 7 and 9, see Figs. 18-20. All of them are proportional to the mass of the charged leptons squared and with two derivatives. Naively, one could expect that the dominant contribution comes from the dimension 5 diagrams. However, as we are considering 4π couplings and lowering the new physics scale as much as possible, all 12 diagrams could be in principle relevant.
We shall show in detail how we derived the integral of the first diagram in Fig. 20 as an example, denoted as I (9) 1 , and give just the results for the rest. We chose this diagram as it gives a similar prefactor as in the original work [10]. After electroweak symmetry breaking (EWSB), we rotate the diagram to the mass basis, see Fig. 21. H + is the Goldstone boson associated to W µ , which appears explicitly with mass m W in the Feynman-'t Hooft gauge. H + are the two mass eigenstates with eigenvalues m 2 + coming form the mixing of S + and η + , which can be trivially diagonalized by a 2 × 2 rotation matrix R H + with angle β. η R,I are the CP-even and CP-odd components of η 0 , with masses m 2 R,I = M 2 η ∓ 1 2 λ 5 v 2 . Defining k ≡ (16π 2 ) d 4 k/(2π) 3 and assigning momenta in the loop, the integral of the diagram in Fig. 21 in the mass insertion approximation is given by  Figure 19: Dimension 7 mass diagrams in the gauge basis. W µ couples with a derivative to the scalars, so every diagram has the same number of derivatives. We will use the notation I where a = R, I and we have neglected the masses of the charged SM fermions. The sum over free indices can be explicitly done enlarging the denominator of the integral. For example, Defining ∆m 2 0 = m 2 R − m 2 I and ∆m 2 + = m 2 +1 − m 2 +2 , Eq. (49) can be written as, with I  with m ρ ++ , Finally, including the corresponding couplings from the potential Eq. (17), the expression for the diagram in Fig. 21 reads where the Yukawa h and the SM charged fermion masses have been omitted. The computation of the rest of the diagrams in Figs. 18-20 is very similar to the example shown. We only give the results here and omit their calculation. The function F Cocktail in Eq. (19) is given by the sum of the different contributions from the diagrams, i.e.
For simplicity we introduced the notation where each numerator, associated to the derivatives depicted in the diagrams, is defined as while for the denominators, with the mass ratios x defined in Eq. (52). The integrals in Eq. (67) are evaluated numerically using pySecDec.
We now discuss the maximization of F Cocktail . We are interested in this case, because the Yukawa h in Eq. (22) is inversely proportional to F Cocktail , and we want to explore the parameter space where h is small enough to be perturbative and avoid CLFV constraints. In general, every integral I (a) i gets larger for smaller masses or, equivalently, smaller ratios. We shall fix a limit of 100 GeV on the scalar masses of η R,I and H + , and 800 GeV for the doubly charged singlet ρ ++ , see Sec. 3 for details. We set the dimensionless couplings λ 5 , λ ηH , and κ to 4π. For the dimensionful µ couplings we impose the limits µ 1 < 4 max[m +1 , m +2 ] and µ 2 < 4 max[m +1 , m +2 , m ρ ++ ], required to avoid the radiative generation of negative quartic scalar couplings [51].

A.2 KNT model
The   Fig. 6. Neglecting the SM charged fermion masses, one finds which is a dimensionless function of the ratios, Note that Eq. (70) is simple enough to be easily decomposed in terms of 3-loop master integrals [70]. Due to the repetitions of the momenta in the denominator, using relation (22) from [8], one has .
As the second and third terms are identical under the exchange of k 1 and k 2 , one can finally write F KNT in terms of a combination of the master integral G integral given in [70]. The resulting expression is The integral has an analytical expression for x 1i = x 2i = 1.
About the maximum value of F KNT , we proceeded analogously to the cocktail model. In this case, we maximized F KNT /M N i , since the neutrino mass matrix in Eq. (29) is proportional to this ratio. We set a lower limit on the mass of the singly charged scalars of 100 GeV and let M N i = M N free. We found that the maximum is around F KNT 60 with m S 1 = m S 2 = 100 GeV and M N i = 840 GeV. Figure 23: Neutrino mass diagrams for the AKS model in the gauge basis.

A.3 AKS model
In this case, there exists two non-equivalent diagrams shown in Fig. 23, which differ by the crossing of the internal S-lines. F AKS is then the sum of the integrals from both diagrams with the correct normalization, given in Eq. (45),