Charged lepton flavor violating processes in the Grimus-Neufeld model

Charged Lepton Flavour Violating (cLFV) decays constrain the relationship between the neutrino and the scalar sectors of the Grimus-Neufeld model (GNM), an appealing minimal model of neutrino masses. It turns out, that in the scenario, where the seesaw scale is lower than the electroweak one, cLFV is completely defined by the new Yukawa interactions between the additional single heavy Majorana neutrino, the second Higgs doublet and the lepton doublets. Therefore, we derive a useful parameterization for the Yukawa couplings which reproduces by construction the correct PMNS matrix and the correct neutrino masses for both Normal and Inverted ordering at one-loop level. We embed this scenario in the $\texttt{FlexibleSUSY}$ spectrum-generator generator to perform parameter scans. Focusing on the tiny seesaw scale, we show that current $\mu\to e\gamma$ limits provide significant constraints on the scalar sector, and we evaluate the impact of future cLFV $\tau$-decay searches for the cases of discovery or non-discovery. The tiny seesaw scale makes the neutrino sector and the cLFV processes in the GNM similar to the scotogenic and the scoto-seesaw models, so we provide constraints for these models as well.


Introduction
Many years have passed since the discovery of neutrino oscillations [1][2][3][4], yet massive neutrinos are still not in the Standard model (SM). That is not surprising: it is extremely hard to either confirm or exclude all the possible mechanisms that generate neutrino masses due to their weak impact on the sectors that we can actually see in the experiments, see e.g. studies of the low energy effects of the seesaw mechanisms in [5,6]. However, having in mind the impressive precision of current Charged Lepton Flavour Violation (cLFV) experiments [7][8][9][10] some of the scenarios can actually lead to possible signatures and/or restrictions on the parameter spaces, and will become even more restricting in the future, given planned experiments [11][12][13][14][15]. Hence, it makes sense to look at the most constraining scenarios to narrow down the possible space of the neutrino mass mechanisms. The simplest neutrino mass mechanisms can be put into three categories: inducing neutrino masses at tree-level (e.g. all types of seesaw [16][17][18]), generating them at loops (see e.g. [19][20][21]), and combining these mechanisms (e.g. [22,23]). These mass mechanisms can be embedded into more exotic models, for example, the seesaw mechanism can arise from extra dimensions [24,25], or radiative neutrino generation can be realized with composite Higgses [26,27]. While virtually all scenarios can fit the current experimental constraints, it is far harder to claim the strictly excluded regions unambiguously, especially when one studies more general scenarios with an overwhelming number of free parameters (see e.g. [28]). A look at less general models, with more exposed and isolated effects on the cLFV from the neutrino sector, naturally provides a solution to this problem.
The scotogenic model [21] is probably the most popular one in this respect, partially due to a manageable amount of free parameters. With an imposed Z 2 symmetry, it has the scalar sector of the inert doublet model (IDM) [29]. The cLFV in the scotogenic model comes purely from the radiative contributions of heavy neutrinos and scalar dark matter candidates [30,31]. Thus, the neutrino sector acts as a bridge between the dark scalar sector and cLFV.
There is an interesting variation, called the scoto-seesaw model [32,33], with even fewer degrees of freedom introduced. In both models, scotogenic and scoto-seesaw, one obtains larger contributions to the cLFV for lower heavy neutrino masses. This potentially implies stronger restrictions on the scalar sector. We define a tiny seesaw scale to be lower than the electroweak scale (studied in e.g. [34][35][36][37][38]), to discriminate between the low seesaw scale that is sometimes considered to include the TeV scale (e.g. [39][40][41]). Then, in the tiny seesaw scale, the scotogenic and the scoto-seesaw models give essentially the same contribution to cLFV as we get in the Grimus-Neufeld model (GNM) [22] -a seesaw extended two-Higgsdoublet model (2HDM). This connection is caused by an approximate Z 2 symmetry in the Yukawa sector [42] (discussed in the next section), which makes the radiative neutrino mass generation in the GNM similar to the one in the scotogenic and scoto-seesaw models.
Global symmetries, such as Z 2 , are extremely useful for classifying the scenarios of beyond Standard model (BSM) physics. Yet they are not something fundamental. If one imposes a global symmetry only on one part of the Lagrangian, while breaking at some other part of it, at some loop order it eventually leads to breakage of the imposed symmetry in the sector of interest too. From one point of view, one might conclude that imposing a global symmetry is only consistent if it is done on the full Lagrangian. On the other hand, if one looks at global symmetries as something that just helps us to classify possible scenarios of BSM physics, there is no real need to insist on satisfying it exactly at all loop levels. For example, the CP-symmetric 2HDM potential might get CP-violating corrections from the quark sector, as discussed in [43]. There are many studies on CP-conserving 2HDM potential, which nevertheless are meaningful and explore an important possible scenario of the scalar sector. This naturally leads to a question: what if some parameters in the Lagrangian are so small, that we can almost see the global symmetry? This generalizes the concept of the global symmetry of the Lagrangian by allowing small deviations from zero of the symmetry breaking terms. Following this logic, we will loosely define the approximate symmetry to be the case in which the Lagrangian parameters that explicitly break the global symmetry are not necessarily zero, but are very small (in our case, smaller than O(10 −7 ) ). The GNM satisfies this requirement in the tiny seesaw limit.
By itself, the GNM is an appealing model of neutrino masses which postulates the existence of only a single heavy neutrino to accommodate neutrino masses and mixings, while the other models require at least two additional fields. The GNM is then particularly attractive in the tiny seesaw parameter region: it has less particles and, because of the approximate (but not exact) Z 2 symmetry, it can express the main phenomenological features of both the scotogenic and the scoto-seesaw models of the cLFV processes and the neutrino sector. Investigating this parameter region we will address the following question in this paper: Which restrictions are imposed by the cLFV on the scalar sector in the GNM?
In section 2, we give an overview of the GNM by introducing the scalar and Yukawa sectors. Then we present the one-loop neutrino mass calculation. It leads to a parameterization of the flavor basis Yukawa couplings, similar to the Casas-Ibarra parameterization [44] that automatically reproduces neutrino masses and mixings at one-loop level. We close the section by noting the existing special parameter points in this parameterization, the checks, the numerical stability, and the limitations of our study, using FlexibleSUSY [45][46][47]. In section 3 we present the analytical expressions for cLFV processes. In section 4 we give a short recap of the most important parameters from the Yukawa sector and the scalar sector that control the branching ratios of cLFV processes and we lay out the strategy of our phenomenological study. We then pursue this strategy in section 5, give an interpretation of the results, and conclude our study in section 6. Technical details, such as the explicit derivation of the parametrization, peculiarities of the minimal free parameter set in the Yukawa sector, and the numerical values of some special cases can be found in the appendices.

The Grimus-Neufeld model
The GNM is a general 2HDM extended with one single sterile neutrino. In the limit of the tiny seesaw scale discussed below, the Yukawa sector is approximately Z 2 symmetric [42], thus its predictions for cLFV become similar to the scotogenic and scoto-seesaw models. Scotogenic and the scoto-seesaw models also predict scalar Dark matter (DM), which is a direct consequence of the Z 2 symmetry. The approximate Z 2 symmetry in the GNM will make the potential scalar DM candidate in the GNM not as stable. Yet to be sure if GNM has a viable DM candidate in this scenario, a dedicated analysis is needed, which is beyond the scope of this work. We do not discuss DM in this paper.

Scalar and Yukawa sectors
In principle, cLFV in the GNM can be analyzed for a general 2HDM scalar sector. However, to highlight the similarities between models we consider the scalar potential to be Z 2 symmetric. The Higgs sector of the model contains two Higgs doublets H 1,2 . The potential takes the form: Only the first Higgs doublet H 1 acquires a vacuum expectation value (VEV) v, and we parametrize the two doublets as: .
All the SM particles, including H 1 , are assigned an even parity under the Z 2 symmetry, while H 2 and the additional sterile neutrino N are odd. 1 The latter enters the Lagrangian with, in general, complex Majorana mass term M and new Yukawa-like coupling Y (i) j to the SU (2) lepton doublets ℓ j , where j denotes the generation: 3) The matrix ǫ = iσ 2 combines the two doublets to an SU (2) invariant product. The Majorana mass M is made real by adjusting the phase of N . Terms with the neutrino Yukawa couplings Y (1) to the first Higgs doublet in eq. (2.3) explicitly break the Z 2 symmetry. The SM-like Z 2 -preserving Yukawa sector reads as: Note the usage of H 1 = ǫH * 1 allowing for another SU (2) invariant product with opposite electric charge. In the flavor basis charged lepton masses are defined by Y f f = √ 2m f /v, f = e, µ, τ , which is also the mass eigenstate basis for charged leptons, but not the mass eigenstate basis for neutrinos.
When the Higgs acquires a VEV, the Lagrangian in eq. (2.3) leads to the two nonvanishing Majorana masses for neutrinos, light m 3 and heavy m 4 (in mass eigenstate basis): where which is the usual type-I seesaw mechanism relating the Yukawa coupling Y (1) i , the light neutrino mass m 3 , and the seesaw scale m 4 ≫ m 3 . The m 3 mass is associated with a light neutrino mass and is of O(0.01 eV). If y → 0 for fixed m 4 then m 3 → 0, and the seesaw mechanism is "turned off" for the exact Z 2 symmetry. For y = 0, the seesaw mechanism applies, and for fixed m 3 , eq. (2.6) can be treated as a relationship between the seesaw scale and the Z 2 breaking parameter, allowing us to define an approximate Z 2 symmetry in the GNM by the inequalities: (2.7) The restriction of eq. (2.7) ensures that the Z 2 -breaking Yukawa coupling is at least an order of magnitude smaller than the Yukawa coupling of the electron.
There are at least two non-vanishing masses for light neutrinos. In the GNM the mass of the other SM neutrino is generated at one-loop level via interactions with H 2 . Motivated by a broken Peccei-Quinn symmetry, small λ 5 values successfully generate radiative neutrino masses with relatively large Yukawa couplings to the second Higgs doublet, namely Y (2) . Just as the limit y → 0 turns off the seesaw mechanism, setting λ 5 → 0 turns off the radiative mass generation leaving one of the needed neutrino masses unexplained. In the GNM the lightest neutrino is massless.
The Z 2 symmetry in the Yukawa sector leads to the similarity between the GNM, the scotogenic, and the scoto-seesaw models in the cLFV phenomenology, even though the Yukawa sectors are slightly different. The scotogenic model [21] has an exact Z 2 symmetry, forcing y → 0 and turning off the seesaw mechanism at the cost of adding 2 additional heavy, Z 2 -odd neutrino states. Then the Yukawa couplings Y (2) , which lead to radiative mass generation, are contained in a 3 × 3 matrix instead of a 3 × 1 as in the GNM.
In the scoto-seesaw model [32,33], one has an exact Z 2 symmetry turning off the seesaw mechanism for the Z 2 -odd N , at the cost of adding one Z 2 -even sterile neutrino, for which, only the seesaw mechanism is allowed. This effectively gives one additional independent parameter, the mass of the Z 2 -even neutrino in contrast to the GNM. In turn, this allows to control the sizes of parameters that enter the radiative and the seesaw mass mechanisms independently from each other in the scoto-seesaw model, while they are related in GNM.
In the GNM, the two SM-like massive neutrino mass states mix in general. To get a convenient parameterization for the Yukawa couplings, one solves the equations for neutrino masses and mixings at one-loop level. We further describe the needed rotations, convenient basis and neutrino mass generation in GNM leading to this parameterization in detail in the following sub-sections.

Rotations for neutrinos
The GNM contains four neutrino states, comprised of the three neutrino components ν i of the lepton doublets, and the single sterile neutrino N . At tree level, neutrino masses arise from the Majorana mass term M and the Yukawa coupling vector Y (1) , coupling the neutrino states to the VEV. These tree-level mass terms of eq.  Figure 1: A sequence of rotations for neutrino mass matrix M F ν . Explicit entries show quantities used in the paper. We denote an expression that is zero exactly at tree-level by 0 0ℓ . With 0 1ℓ we mean an approximated one-loop zero, where terms proportional to y 2 × loop and m 3 m 4 × loop are neglected. The matrixΣ contains non-vanishing expressions for the processes ν ′ i → ν ′ j at one-loop level, which are not further suppressed by y 2 or m 3 m 4 . This decomposition is useful because each rotation leads to a separate physical consequence or highlights important details about the Yukawa couplings. The sequence of these three rotations and their effect on the mass matrix is schematically shown in figure 1. In the following, we will discuss the diagonalization steps in more detail. It will be useful to introduce the decomposition of the relevant 4 × 4 matrices as: where U , V and R are 3× 3 unitary matrices andR is a 2× 2 unitary matrix. We will define the matrixS in eq. (2.12) and explain the approximate equality right after eq. (2.14).
In the Lagrangian of eq. (2.3) new complex Yukawas Y (i) are introduced, which carry 12 degrees of freedom in total. As the first step of the diagonalization process we use the freedom to rotate flavor neutrinos ν α to choose a convenient basis. This is done by the 2 Our conventions for unitary rotations are the ones of ref. [48], as implemented in FlexibleSUSY: Figure 2: Diagrams, contributing to theΣ mass matrix. Arrows show the flow of chirality.
matrix V in the following way: This phase, α 1 , is used to absorb one of the Majorana phases η i of the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix. Additionally, 3 degrees of freedom in V can be absorbed into the phase definitions of neutrinos, leading to 5 physical degrees of freedom in V , which are later related to the 3 angles, 1 CP phase, and 1 Majorana phase of the PMNS matrix. 2 of the 4 degrees of freedom of the rhs. of eq. (2.10) will be related to neutrino masses. It is clear from eq. (2.10) that one of the neutrinos doesn't interact with the Higgs bosons after the rotationṼ , which leads to vanishing contributions to its mass at tree level and at one-loop level. These vanishing elements of the mass matrix are shown by 0 1ℓ in the top and left outer entries of the matrices in figure 1.
At the second step of figure 1, the tree-level mass matrix has a typical seesaw structure. Correspondingly, the second diagonalization matrixS is a tree-level seesaw transformation, which yields a diagonal mass matrix at tree level with two non-vanishing tree-level mass eigenvalues m 3 and m 4 , related by eq. (2.6). The appropriate seesaw rotationS can be written as:S After this diagonalization step, we take into account one-loop corrections to the masses. Since one of the neutrinos does not couple to the Higgs bosons, the one-loop neutrino mass corrections lead to the block structure indicated in the third matrix of figure 1: the two non-vanishing blocks at the level of our approximation are the 2 × 2 mass matrix block Σ, which contains the tree-level eigenvalue m 3 and one-loop self-energy corrections, and the heavy seesaw neutrino mass m 4 . All other entries are either generated at higher loops or are neglected in one-loop diagrams as being either proportional to y 2 × loop, which is extremely small in our case (see eq. (2.7)), or to m 3 m 4 × loop, which has an additional seesaw suppression factor. As will be shown in the next section, the last rotationR diagonalizesΣ and can be used to parameterize the Yukawas couplings Y (i′′) in the mass eigenstate basis. Its form will be described in eq. (2.21).
There are in total three non-zero masses for neutrinos at one-loop level. We write them in the ascending order as:m where we consider the loop corrections for the heaviest mass to be negligible. Also, m 4 ≫ m 3 leads toŨ =RSṼ ≈RṼ ⇒ U ≈ RV, (2.14) which means that neutrino masses and mixings can be related to the so-called 3ν mixing paradigm |Ũ i4 | ≪ 1, which was indicated in the definition ofŨ in eq. (2.9). The PMNS matrix is defined as [49]: where η 1 and η 2 are unknown Majorana phases, s ij = sin θ ij , and c ij = cos θ ij , see eq. (2.18).
In the GNM, the lightest neutrino has a vanishing mass and the phase η 1 can be absorbed into a redefinition of the corresponding field, as shown in eq. (2.11). Hence only η 2 is physical. To keep our study simple, and since there is no evidence of η 2 in the PMNS matrix so far, we set it to zero. The zero-mass lightest neutrino also implies a lower bound on neutrinoless double-beta decay [50], which is however left out of the scope of this paper. The pole masses in eq. (2.13) are different for Normal ordering (NO) and Inverted ordering (IO). Since the mass of the lightest neutrino vanishes, the measured mass squared differences from the neutrino oscillation experiments determine the actual neutrino masses. Therefore, eq. (2.8) connects the PMNS matrix with U , which we can summarize as: where the sign convention for ∆m 2 32 is explicitly avoided, and the mass ordering of IO is related to the usual 3ν conventions by: We take as numerical values in our code from ref. [49]: where δ CP is taken from the intersection of measured 1σ-regions for NO and IO scenarios.

Neutrino mass matrix
The Yukawa interactions introduced in eq. (2.3) can be parameterized by four real parameters, as shown in eq. (2.10). One can connect them with the definition of the unitary 2 × 2 rotationR and use the compact parameters of the latter instead. This section provides the required relations. The light neutrino mass matrixΣ disentangles from the heavy one due to eq. (2.14), see ref. [51] and figure 1. The interactions with the doublet H 1 are proportional to the Yukawa coupling y ∼ |Y (1) | in eq. (2.7), which we assume to be vanishingly small due to the approximate Z 2 symmetry. The relevant contributions toΣ are shown in figure 2: with the loop contribution similar to the scotogenic model [21], and the definition of Λ as in ref. [52]. 3 The diagonalization ofΣ is done via the Takagi decomposition: which restricts the matrixR, as shown in appendix A, due to the zero determinant of the one-loop correction term. For the rotation itself Murnaghan's parameterization is used: (2.21) The following ranges of angles and phases define the rotation in a unique way: The parameter ranges that uniquely describe cLFV ratios are smaller: The four degrees of freedomm 3 , d and complex d ′ -are replaced by the two neutrino pole masses m pole 2,3 and the rotation matrix parameters r and ω 22 . The other parameters of R are fixed by eq. (2.20). In this paper we take the point of view that very large one-loop corrections for m 3 are fine tuned and unnatural: hence we only analyze parameter regions where they do not exceed 50%, which corresponds to the following range for m 3 : (

2.24)
It is convenient to define the abbreviation t 32 , representing the ratio between the two physical neutrino masses: The expression for Λ given in ref. [30] is two times larger, which looks like a direct equivalence to eq. (11) of ref. [21]. However, the definitions of the Yukawa couplings in these references differ by √ 2, which should lead to additional factor of one half for Λ in ref. [30]. Hence, there might be a typo in the neutrino Yukawa couplings in ref. [30]: they are by a factor of √ 2 smaller then required.  .27), except for the red striped region, where |z| < 1 2 , which we exclude as being fine-tuned. These red striped areas will also be colored white in later plots, as we will not consider this parameter region for the discussed exclusion criteria. The points where the second Yukawa coupling to electrons, muons, or taus vanish exactly are also shown for both hierarchies.
It is also convenient to define the parameter z by z(r , ω 22 , ω 32 ) = cos 2 r e 2iω 22 + t 32 sin 2 r e 2iω 32 . (2.26) Using (A.5) we see that z corresponds to the relative loop contribution to the neutrino mass m 3 , and we can express eq. (2.24) in terms of these parameters:  .27) for ω 32 to be used in the parameterization. However, these solutions lead to physically equivalent Yukawa couplings, as shown in appendix B, thus we use The region for the parameters r and ω 22 satisfying the constraint of eq. (2.27) is shown in figure 3 as colored regions, while white regions are disallowed.
With these definitions, we can determine the Yukawa couplings Y (i′′) in the one-loop mass eigenstate basis, which reproduce the PMNS matrix and the neutrino masses by construction. The following Yukawa couplings follow from eq. (2.20): Changing the sign of Λ changes the phase of Y (2′′) . But this phase cancels in the calculation of cLFV decays. This means that cLFV observables restrict only the absolute value of Λ. The relation between mass eigenstate Yukawa couplings Y (i′′) of eqs. (2.30, 2.29) and the flavor ones Y (i) is given by: The rotation U is related to the PMNS matrix in eq. (2.16). With the Yukawa coupling of eq. (2.29) this resembles the Casas-Ibarra parameterization, ref. [44], adapted for the scotogenic model, ref. [30]. Setting m pole 3 → 0 in eq. (2.29) realizes the exact equality.

Vanishing Y (2)
Neutrinos contribute to cLFV processes via Yukawa couplings to the second Higgs doublet, Y (2) . Hence, the vanishing Yukawa coupling leads to a vanishing neutrino contribution to the process of interest. It turns out, there are solutions that simultaneously give a vanishing Yukawa coupling in the flavor basis and reproduce neutrino masses and mixings. Using eqs. (2.21, 2.29, 2.31) and setting Y (2) f = 0 for a given flavor f one gets: The solutions of eq. (2.32) are not backed by any fundamental reasoning, are not stable and are fine-tuned. Nevertheless, these solutions do exist and, if the small parameter ranges around them are included, they will lead to weaker constraints on the scalar sector for cLFV. Thus, if we discard some tiny parameter region around these points, we would make our main result, a constraint on the scalar sector from cLFV, dependent on our own choice of the size of the region. Instead, we choose a more agnostic approach: we include it, discuss it in more detail, and provide a suggestion of how one can define a "more likely" scenario away from these regions and a "fine-tuned" one around these regions. In this way, we completely cover the full parameter space of the model and do not make our main result dependent on additional assumptions.

Reproducing neutrino masses and mixings
To realize the numerical scans over model parameters we used the FlexibleSUSY [45][46][47] spectrum-generator generator, including the extension NPointFunctions [53]. It allowed us to straightforwardly implement the non-trivial mass-generation mechanism for neutrino masses in the GNM, combining tree-and loop-level generated masses, as well as the experimentally measured PMNS mixing matrix. In particular, we created SARAH [54,55] model files for the realization of the GNM studied in this paper. In addition, we designed the FlexibleSUSY model file incorporating the parameterization of the parameter space developed in section 2. 4 This working setup also resulted in an independent cross-check of the consistency of our analytical one-loop parameterization of the neutrino sector.
However, using the current implementation of the code FlexibleSUSY we spotted a bug in the neutrino pole mass calculation, where the couplings in self-energies taken from SARAH were conjugated, compared to analytical expressions both done by hand and with FeynArts [56] and FormCalc [57]. As a workaround, the conjugation of self-energies in the neutrino pole mass calculation in FlexibleSUSY was applied.
For the studied tiny seesaw parameter region, the elements of the PMNS matrix and the output neutrino masses at one-loop are consistent within 1% accuracy. However, for some limiting cases that lie beyond the region of our interest, the description ceases to be accurate enough.
Since we study a tiny seesaw region, scale differences in the neutrino sector are limited only to 12 orders of magnitude (m 4 ≈ 10 GeV vs. m 2 ≈ 10 −11 GeV). This helps us in numerical stability. In fact, we checked that for a seesaw scale of m 4 > 10 GeV, the neutrino mixing matrix in FlexibleSUSY becomes inaccurate (>1%), while for the neutrino masses the scale is higher m 4 > 100 GeV. We find numerically stable and correct neutrino masses and mixings, allowing up to 1% deviation in the output of FlexibleSUSY, for Λ > m pole 3 ≈ 5 · 10 −11 GeV . (2.33) It is interesting to note that this limiting value gives the largest Yukawa coupling with value Y (2) ≈ O(1). Hence, the Yukawa values higher than O(1) and up to a perturbativity limit do not accurately reproduce the neutrino spectrum in FlexibleSUSY for this model. For our study, however, we will not reach this scenario since cLFV gives stronger constraints in general. This means that all our results are consistent, cross-checked, and reproducible with FlexibleSUSY. 5

Charged Lepton Flavor Violating processes
The new Yukawa interactions are not only responsible for the generation of neutrino masses, but they also give rise to cLFV. In this section we provide the theoretical formulas for the amplitudes of two-body decays l i → l j γ, µ → e conversion (valid for all nuclei), and three-body decays l i → l j l k l c k . First, we present Feynman diagrams and the amplitudes, specifying relevant and negligible ones. Later, the amplitudes are combined into the decay (and conversion) rates. The numerical calculations are also implemented in FlexibleSUSY, using the additional extension NPointFunctions [53].

Coefficients
Let us start from the simplest case of the penguin contributions, see the Feynman diagrams in figure 4. We express the amplitude of the flavor-changing decay l i → l j γ into an off-shell photon with outgoing momenta q = p i − p j as in ref. [58]: For the considered scenario in the GNM, the dominant photon contribution is given by one-loop diagrams of figures 4b-4d with virtual charged Higgs boson H − and virtual Majorana neutrino ν ′′ 4 exchange. In these diagrams, the lepton flavor transition is mediated by the new Yukawa coupling Y (2) .
Other contributions are negligible due to the following reasons. The impact of Goldstone bosons G − W is suppressed by the Z 2 -breaking Yukawa coupling |Y (1) | ≪ 1. The impact of vector bosons is suppressed due to the Glashow-Iliopoulos-Maiani (GIM) mechanism which leads to a m 2 ν ′′ i /m 2 W factor and to terms proportional to |Ũ i4 | ≪ 1. Also, in our scenario, for other cLFV processes under consideration, one can always neglect the electron Yukawa coupling Y (1) due to its small magnitude relative to Y (2) .
The only non-vanishing coefficients in the decay rate, eq. (3.1), are the following: with the loop functions given in eq. (3.5).
Another class of cLFV diagrams are Z-boson penguins, shown in figure 4. It turns out, that they can be removed from the consideration. The diagram in figure 4a is proportional to the ν ′′ α ν ′′ β Z coupling. Though this coupling does not vanish for the first three generations of neutrinos, flavor-changing couplings to external leptons includeŨ α4 , which leads to a seesaw suppression. The diagrams of figures 4b-4c have non-vanishing couplings but they lead to the zero factor of (B 1 + 2C 00 ) with 3) The last diagram in figure 4d is proportional to Y (1) . Overall, this discussion shows that the impact of Z-penguins is negligible for our scenario. Similarly, the Higgs boson penguin is suppressed by Yukawa couplings to charged leptons. From this follows, that only the photon penguin gives a non-vanishing contribution. Next, we consider box contributions relevant for l i → l j l k l c k , see figures 5a-5d. They sum up to the result:  Figures 4c-4d do not give numerically meaningful contribution, but they have to be included to ensure that UV divergences cancel exactly. All the diagrams that include Z boson are negligible, while the diagram shown in figure 4a is zero for a photon contribution. For the µ → e conversion in presence of a nucleus, there are relevant box diagrams as well, see figures 5e-5f. They lead to a vanishingly small contribution, because in the GNM quarks are coupled to H 1 , leading to a Y (1) suppression in addition to the small Yukawa couplings.
The loop functions used above have the following form (see [59]): (3.5) The relevant limit for us is |x| ≪ 1, and the loop functions are normalized as F i (x → 0) = 1.
Using eqs. (2.29, 2.31) one obtains for m 4 ≪ v m H ± that while the only other parameters that the amplitudes depend on are r and ω 22 . This dependence is only slightly more complicated and can be read out from the Yukawa couplings.
Note, that the single parameter relating the scalar sector to the cLFV two-body decays is |Λ|m 2 H ± . This is also the most important factor for three-body decays because contributions from box diagrams can be neglected for m H ± TeV. In addition, the model predicts a very simple relationship between the two photonic amplitudes: (3.7)

Decay widths
We use the following convention for a covariant derivative: where e > 0 is the charge unit, Q f the electric charge of a corresponding fermion f and A µ the photon field. The total decay width of l i → l j γ simplifies from eq. (3.1) to (see, e.g. [58,59]): The partial decay width with three leptons of the same generation in the final state is (see [60][61][62]) where the u-channel for photon penguins is included via Fierz identities; for boxes, all channels are calculated directly, and the following abbreviations are used: The minus sign for the photon penguin comes from the embedding into a four-fermion amplitude, and the factor 1/2 for the boxes comes from Fierz identities during the matching. For three leptons of different generations in the final state, the expression for the partial decay rate differs [63]: with different matching for boxes leading to an absence of the additional prefactor: For the conversion rate, we use [64] ω µ→e = 4m 5 with dimensionless integrals D and V (p) defined there as well. The minus sign reflects the different definition of the photon field that comes from the comparison of covariant derivatives.

Recap of important parameters
For our study we restrict the parameters: The first inequality rewrites the condition on the seesaw scale of eq. (2.7), the second one leads to charged Higgs boson masses well in reach of the LHC, and also leads to negligible box contributions in almost all of the parameter space and significantly simplifies the discussion. 6 The cLFV ratios are determined by the Yukawa couplings to the second Higgs doublet, Y (2) (ω 22 , r , Λ), which also includes neutrino masses and mixings determined by the oscillation parameters. Since the parameter |Λ| factors out in the amplitudes as in eq. (3.6), all the relevant parameters for cLFV are ω 22 , r, |Λ|m 2 H ± and Λ 2 m 2 H ± . However, the last one comes from the box diagrams, which is out of the scope of this paper due to a typically negligible box contribution impact for cLFV with relatively light charged Higgs bosons. Neglecting the box contributions, all the relevant free parameters that enter cLFV then sum up to: The first two are the free parameters that define the Yukawa sector: ω 22 is a phase that is related to CP and Majorana phases at one-loop, while r parameterizes the mixing between seesaw and radiative states. The last parameter of eq. (4.2) relates the scalar and Yukawa sectors. We will refer to it as the photon factor, as it is a factor in front of the amplitudes that include a photon, i.e. in front of A L 1 and A R 2 , as shown in eq. (3.6). The photon factor, |Λ|m 2 H ± , relates the scalar sector to the Yukawa sector via cLFV and it is of most interest. Using eq. (2.19), we write it explicitly: For the tiny seesaw region, |Λ|m 2 H ± can be generalized to the general 2HDM [42], by inclusion of the mixing of the scalar particles in Λ, eq. (2.19). In our study, we give bounds on the photon factor that are independent of the exact form of the scalar potential or its symmetries (Z 2 or CP breaking or not). Note, however, that the Z 2 symmetry breaking Yukawa couplings of charged leptons to the second Higgs doublet in the Higgs basis can in principle alter the cLFV values, if included.
It is instructive to look at the dependence on the parameters of the scalar potential for the sake of physical intuition. The inert scalar potential of eq. (2.1) in the limit of approximately degenerate heavy scalar masses simplifies the photon factor: 7 (4.4) 6 See next section. 7 An exact degeneracy would lead to λ5 → 0 and vanishing radiative neutrino mass.  Table 1: Current and planned experimental bounds, related to corresponding observables. Data for τ decays for Belle-II was obtained from the figure 189 of ref. [12]. For the µ → e conversion, we consider the Al nucleus.
The non-observation of cLFV then generally leads to a lower bound on the photon factor |Λ|m 2 H ± (and thus λ 5 ) as a function of the seesaw scale. As was mentioned before, λ 5 is a Peccei-Quinn symmetry breaking parameter that is bounded from below by the neutrino sector alone. The cLFV decays then allow us to improve this bound in the tiny seesaw region.
As a last note, we stress again that the cLFV rates are related to the scalar sector only via a single parameter, |Λ|m 2 H ± . Hence, our results do not depend on how exactly all the parameters in the potential are realized. While we used the scalar potential of the Inert Doublet model (IDM), motivated by an approximate Z 2 symmetry, and assumed that the values of scalar masses are close to each other to get the specific interpretation of a |Λ|m 2 H ± in eq. (4.4), one is, in general, free to use any scenario of the 2HDM potential, as already mentioned before. To construct a consistent potential, one can look at e.g. [65][66][67][68][69][70] for detailed studies of 2HDM or at [71] for constraints in IDM specifically.

Discussion of the relative importance of different decay modes
All the cLFV experiments, which potentially can constrain the GNM, are shown in table 1. Below, we discuss the importance of all these decay modes to single out the most constraining experiments for the GNM.
Typically, the branching ratios for three-body decays are dominated by photonic contributions, while others -boxes, Z and Higgs penguins -are negligible in the parameter region defined by eq. (4.1). The tiny seesaw scale, defined as m 4 ≪ v m H ± , further leads to a fixed ratio of A L 1 ≈ 2/3A R 2 . We call this regime photon dominance (in contrast to dipole dominance, which assumes A L 1 ≪ A R 2 ): where α = e 2 /(4π). The first term of eq. (5.1) is a correction to the dipole dominance. This correction amounts to O(10%) for µ → 3e. The box contribution can be increased for lower |Λ| values, as follows from eq. (3.6). However, two-body decays constrain the |Λ|m 2 H ± factor from below. Taking this minimum value of |Λ|m 2 H ± , a lower value for |Λ| translates into a higher value for m H ± . We concentrate on the lighter masses of the charged Higgs in the paper, see eq. (4.1), where deviations from eq. (5.1) are small.
We find that the box contributions are generally negligible in all of the ω 22 − r plane for µ → 3e, except for the very small region, where the rate for µ → eγ drops sharply close to the point where Y (2) µ = 0. This is because the constraint on |Λ|m 2 H ± from µ → eγ is much looser around that parameter point than in all the rest of the parameter region. Since it is contained in a tiny enough area in the ω 22 − r plane and both, BR(µ → 3e) and BR(µ → eγ) go to zero at Y (2) µ = 0, it gives negligible modifications to the information that the two-body decays provide.
The only parameter region where τ → 3µ can be expected to be observed in Belle-II is around Y (2) e = 0; such an observation can happen only if τ → µγ is also seen. The value of τ → 3µ can in fact deviate more significantly from photon dominance, but we find that the current and planned experimental sensitivities still leave this process phenomenologically irrelevant. Other three-body processes are not expected to be seen in Belle-II at all. As a result of these checks, we conclude that in the parameter region of our study, the three-body decays are of minor importance compared to the two-body ones.
The second phases of experiments, Mu3e-II and COMET-II, will enhance the importance of the corresponding processes. Their branching ratios will be fixed by the photon dominance contributions in the studied parameter region as box contributions can be neglected. In the present paper, we do not consider these longer-term improvements. For the next achievable improvement of µ → e conversion at COMET [14] the following relation holds: which makes it less restrictive than µ → eγ. From now on, we will concentrate on the study of two-body decays as the most constraining decay modes for the GNM. Note, however, that two-body decays can vanish if the corresponding Yukawa couplings vanish. These points in the parameter space do exist, as described in section 2.4. Nevertheless, as can be seen from figure 3, there is no such point in the ω 22 − r plane in which two of the Yukawa couplings in flavor basis vanish simultaneously. This means that all three two-body decay experiments have to be combined to give a strict bound on |Λ|m 2 H ± and, hence, we further study all three of them.

Current and planned restrictions on |Λ|m 2 H ±
The lower bound for |Λ|m 2 H ± from non-observation of cLFV depends on the parameters r and ω 22 . The current lower bounds from µ → eγ on |Λ|m 2 H ± for NO and IO are given as contour plots in figure 6. The most important observable is µ → eγ. It gives the tightest bounds almost everywhere, except for the sharp minima around the points, where the   corresponding Yukawa couplings vanish. In those small regions, either τ → eγ or τ → µγ constrains the photon factor |Λ|m 2 H ± , which results, for the current experimental limits (first column of table 1), in the lower bounds on the photon factor, shown in table 2 (also see figure 8 for current and planned future bounds).
outcomes of the experiments: where superscripts "current" and "planned" are again used for the reach-in branching ratio of the current or planned experiments. The corresponding values are given in the first and The parameter regions where the observation of τ → eγ and τ → µγ in the ω 22 − r plane is possible are shown in figure 7. They also provide the maximal allowed photon factor by the values given in table 3. The lower bounds from table 2 together with the upper ones  from table 3 specify the allowed range for the photon factor |Λ|m 2 H ± in case of these observations, while figure 7 shows the allowed regions in the ω 22 − r plane: the dashed areas. These dashed areas are disjoint, meaning that the GNM predicts, that τ → eγ and τ → µγ cannot be seen together. If µ → eγ is observed, then figure 6, scaled with eq. (5.3) for the observed branching ratio, gives the values of |Λ|m 2 H ± in the ω 22 − r plane.

Nothing is observed: future absolute bounds on |Λ|m 2
H ± To generalize our results and to give a more definite answer of how the combination of the neutrino sector and cLFV constrains the scalar sector, we look at the absolute bound on the photon factor, |Λ|m 2 H ± . By absolute bound we mean the minimum value of |Λ|m 2 H ± such that the whole ω 22 − r plane, consistent with neutrino masses and mixings, is excluded by cLFV observables. The absolute limits for current experimental bounds are shown in the first row of table 2 and are determined by τ → eγ experiments at Y In general, µ → eγ is the strongest experiment in most of the ω 22 − r plane in the nearest future. However, the allowed minimal |Λ|m 2 H ± is mainly determined by the weaker bounds on τ → eγ or τ → µγ. This is because one always has two points in the parameter space where Y (2) e,µ = 0 and, hence, the branching ratio for µ → eγ vanishes. In order to further improve the absolute bound on |Λ|m 2 H ± one therefore has to improve the experimental sensitivities on τ → eγ and τ → µγ. The combination of branching ratios of these exper-iments, required to completely exclude the ω 22 − r plane for specific values of |Λ|m 2 H ± , is shown in figure 8.
The yellow borders specify precise values of branching ratios that have to be applied, i.e. if the future τ → eγ and τ → µγ sensitivities are inside one particular yellow rectangle, the absolute lower bound on |Λ|m 2 H ± is at least as good as indicated in the legend. The yellow borders are obtained without approximations and are determined by varying the branching ratio of the τ decays over the regions in the ω 22 − r plane that are not excluded by µ → eγ. The blue rectangles are similar, but they show the hypothetical case in which µ → eγ is of infinite precision, and hence only the points with BR(µ → eγ) = 0 are allowed by µ → eγ. All future improvements in µ → eγ experiments for a specific value of |Λ|m 2 H ± reside in the area between the closest yellow and blue rectangles. One can see, that the size of the area between the rectangles is small (hence for this analysis improvements in µ → eγ are of minor importance) but grows with the values of the photon factor. This behavior continues to higher values of |Λ|m 2 H ± than shown in the plot, up to the moment when µ → eγ allows the point where Y (2) τ = 0, as shown as critical case in figure 6. To exclude even larger |Λ|m 2 H ± values, the improvement of bounds for µ → eγ is required. For a different µ → eγ precision, figure 8 can be used with the following rescaling of the photon factor and axes: However, a high precision of τ experiments is required so that this case lies out of the bounds of figure 8. Taking into account the current and planned bounds for τ → eγ and τ → µγ, one concludes that only τ → eγ improvements are important for the absolute bound in the near future.

Interpretation of results: example of limits on λ 5 and relations to scotoseesaw and scotogenic models
The absolute bounds on the photon factor, |Λ|m 2 H ± , correspond to the values of the photon factor, for which the entire ω 22 − r plane is excluded. They are given by τ → eγ in table 2 or figure 8. However, figure 6 shows that the absolute bounds are due to very small areas in the parameter space. These areas are characterized by the vanishing of one of the Yukawa couplings, and only in those regions an observation of τ decays in Belle-II can be expected. We define the "essential part" of the ω 22 − r plane, as the region where the observation of τ decays in Belle-II is not possible i.e. all the plane except the regions shown in colors in figure 7. We then obtain a "typical" bound on the photon factor, i.e. the bound for which the "essential part" of the ω 22 − r plane is excluded. The absolute bound is shown in table 3 and together with this typical bound, defined above, reads as: absolute for NO(IO): |Λ|m 2 H ± > 1.9(4.0) · 10 −6 GeV 3 , typical (no τ → e(or µ)γ expected): |Λ|m 2 H ± 10 −4 GeV 3 . (5.6) To get an intuition of how strong this current bound on the scalar sector actually is, we can go back to the special case of the scalar sector, where eq. (4.4) holds. Then eq. (5.6) translates into bounds on |λ 5 |: absolute for NO(IO): typical (no τ → e(or µ)γ expected): |λ 5 | keV m 4 . (5.7) As we mentioned in the introduction, the GNM in the tiny seesaw scale region is similar to the scoto-seesaw and scotogenic models due to an approximate Z 2 symmetry. To put our results in a more general framework let us consider how our results can be applied to these two models.
The scoto-seesaw model has two heavy neutrinos, one of which is odd under the Z 2 symmetry and the other one is even. One can use our parameterization for the Yukawa couplings, where m 4 in Y (2) is understood as the mass of the Z 2 -odd neutrino, while m 4 in the expression for Y (1) has to be understood as the mass of the Z 2 -even neutrino of the scoto-seesaw model. In case these heavy neutrinos of the scoto-seesaw model are lighter than the scalars, only Y (2) enters the cLFV just as in our case, and hence all our results for cLFV, including all figures, hold for the scoto-seesaw model, too. We also note that our parameterization allows to cover the full parameter region of the scoto-seesaw model, while [32,33] assume only scenarios, where the mixing between the radiative and the seesaw neutrino at one-loop is absent, i.e. r = 0 in our parameterization.
The scotogenic model also has a Z 2 symmetry, yet it does have more parameters. Nevertheless, the Casas-Ibarra parameterization gives the same general behavior of the Yukawa couplings with respect to |Λ|m 2 H ± and hence cLFV ratios like in our case. The difference is that in the scotogenic model |Λ|m 2 H ± is a 3 × 3 matrix and, instead of two parameters, r and ω 22 , one has a general 3 × 3 orthogonal matrix entering the Yukawa couplings. Nevertheless, the same typical behavior for the tiny seesaw scales in the scotogenic model is to be expected. One should also have in mind that in the scotogenic model there could be similar cancellations suppressing the µ → eγ branching ratio like in our case. These cancellations would push the absolute bound on |Λ|m 2 H ± lower. However, from our study we cannot generalize our strict exclusions and predict them for the scotogenic model.
Even though the tiny seesaw scale has never been rigorously studied in the scotogenic model before, we do find the comment about it in [30], which gives us some way to put our results in a more general context. [30] claims that for m 4 ≪ m H ± the value λ 5 = 10 −9 is excluded. We see that the typical bound, given in eq. (5.7), indeed is consistent with this statement for scalar masses of the TeV scale, which are the masses they actually consider. This confirms our expectation that the typical behavior of the models in this parameter region is the same. As said before, our absolute bound from eq. (5.7) cannot be directly applied to the scotogenic model, yet it gives an example of how much special analytical solutions, that can hardly be caught in random parameter scans, can push the strict exclusion limits of the model.

Conclusions
Throughout this paper, we concentrated on the question: how does the neutrino sector of the GNM constrain the scalar sector via cLFV observables? The GNM is a model where neutrino masses are generated with a minimal neutrino sector but a non-minimal Higgs sector, such that new Yukawa couplings can lead to cLFV effects. Guided by the principle "exclude as much as possible", we singled out the scenario where the strongest constraints can be drawn, i.e. when cLFV decays are enhanced. This scenario corresponds to a tiny seesaw scale (lower than the electroweak scale) and a mass of the charged Higgs m H ± slightly above the electroweak scale (up to a TeV).
In this parameter range, box diagrams are negligible, and the study simplifies to the relations between three parameters: two parameters, ω 22 and r, parameterize the new Yukawa sector, see figure 3, and one parameter, the so-called photon factor |Λ|m 2 H ± , includes scalar sector parameters and the seesaw scale, see eq. (4.3).
The non-observation of cLFV then results in lower bounds on the photon factor as a function of ω 22 and r. These lower bounds are shown as contour plots in figure 6 for each neutrino mass ordering. Currently, µ → eγ is the most constraining observable in the largest part of the ω 22 − r plane. There are only two areas in this plane where µ → eγ gives weaker constraints than τ → eγ and τ → µγ, corresponding to areas of approximately zero Yukawa couplings Y In the mentioned special parameter areas, the observation of τ → eγ or τ → µγ in the planned experiments becomes possible and leads to the regions shown in figure 7. The upper limits on |Λ|m 2 H ± for these decays is given in table 3. Hence, if τ → eγ or τ → µγ is observed in the future, the parameter space is drastically reduced, while the observation of both τ → eγ and τ → µγ is excluded in the GNM by the sensitivity of µ → eγ.
To disentangle the discussion of the neutrino sector and the scalar sector, we look at global bounds on |Λ|m 2 H ± , which are independent of the neutrino parameters ω 22 and r. One option is to consider absolute bounds, i.e. bounds for which the complete ω 22 − r plane is strictly excluded. To increase this absolute bound an improvement of the sensitivity in τ → eγ is most important in the near future. This situation is shown in figure 8, where the needed branching ratios of τ → eγ and τ → µγ to exclude the particular value of |Λ|m 2 H ± for the complete parameter space of the model can be extracted. For instance, until the BR(τ → eγ) ≈ 10 −11 (10 −10 ) for NO(IO) is probed, any other improvement of cLFV observables has no effect on the bound on the photon factor |Λ|m 2 H ± . Technically, the absolute bound is driven by the point, where Y (2) µ = 0, see table 2. In the numerical evaluation it is helpful to rely on analytical solutions for such points, provided in section 2.4.
A second option is to consider typical bounds, i.e. bounds on |Λ|m 2 H ± for which the essential part of the ω 22 − r plane is excluded. The absolute and typical bounds are shown in eq. (5.6) and interpreted in terms of the scalar sector parameter λ 5 in eq. (5.7).
All the presented results are directly applicable to the scoto-seesaw model and complement the current studies of [32,33]. Additionally, the presented neutrino mass calculation and the parameterization of the Yukawa couplings are more accurate and do not neglect the mixing between the radiative and the seesaw neutrino states at one loop as in [32,33]. This allows us to cover the parameter space completely.
One cannot apply the absolute bound from the GNM directly to the scotogenic model. In both models, the dependence on the photon factor coincides, but there are more free parameters in the Yukawa sector of the scotogenic model. This leads to unrelated analytical solutions in the limiting case. However, as discussed in section 5.5, the typical bounds of eqs. (5.6, 5.7) are applicable for the typical behavior of the scotogenic model.

A Parametrization for Yukawa couplings
In this appendix we give the step by step derivation of how we can define the Yukawa couplings unambiguously.
We start with defining the loop contribution toΣ, eq. (2.19), as which has determinant zero. For convenience we repeat the definition ofR: we can restrict the components ofR by the condition, that the determinant of the rhs. of eq. (A.1) vanishes: giving us a definition of the phase φ R : W itself, eq. (A.1), produces also additional relations between d, d ′ andR ij : Since Λ is real, see eq. (2.19), we get from eq. (A.6) meaning that z has to be real, and additionally, has to have the same sign as Λ. From eq. (A.5) we get then e 2iφ R = sign(z) = sign(Λ) , (A.8) which allows two values for e iφ R , depending on the sign of Λ: From our limit on the range of m 3 , eq. (2.24), we get also a limit for |z| = m pole 3 /m 3 : (A.10) Using above definitions and conditions, eqs. (A.3, A.5), it turns out that we can rewrite as a tensor product: Motivated by this decomposition we can extend the tensor product definition, using eq. (2.10), andR from eq. (2.9) to express: Analogously, we have to rotate the Yukawa coupling to the first Higgs doublet: using eq. (2.6, A.5).

B Determining the minimal parameter space
In this appendix we derive a minimal region of φ R , r, ω 22 , and ω 32 that covers the whole parameter space to study the unique constraints by cLFV. The dominant contribution of all cLFV processes in our region of interest, that is motivated by the approximate Z 2 symmetry, depends on Y we see the following relation Y (′′) (r ± π) = −Y (′′) (r) ⇒ r ∈ − π 2 , π 2 , (B.1) that reduces the required region of r, as the minus signs are just a different phase, that does not influence the cLFV. In the same way, we can pick in eq. (A.9) the plus sign before the square root, as the different phase does not influence the cLFV.

C Relations between Yukawa couplings of different parameter points
It turns out, that there is a relation between the parameters ofR, that gives the same Yukawa couplings for both solutions, ω 32+ and ω 32− . We also notice that the conditions for the allowed areas in the ω 22 − r plane depend only on the sizes of |r| and |ω 22 |.