Two-body lepton-flavour-violating decays in a 2HDM with soft family-lepton-number breaking

We evaluate the decays $\ell_1^\pm \to \ell_2^\pm \gamma$, $Z \to \ell_1^+ \ell_2^-$, and $h \to \ell_1^+ \ell_2^-$, where $\ell_1$ and $\ell_2$ are charged leptons with different flavours and $h$ is the scalar particle with mass 125.25 GeV, in a two-Higgs-doublet model where all the Yukawa-coupling matrices conserve the lepton flavours but the Majorana mass terms of the right-handed neutrinos break the flavour lepton numbers. We find that (1) the decays $\ell_1^\pm \to \ell_2^\pm \gamma$ require large Yukawa couplings and very light right-handed neutrinos in order to be visible, (2) the decays $Z \to \ell_1^+ \ell_2^-$ will be invisible in all the planned experiments, except in a very restricted range of circumstances, but (3) the decays $h \to \ell_1^+ \ell_2^-$ may be detected in future experiments for rather relaxed sets of input parameters.

1 − 2 , and h → + 1 − 2 , where 1 and 2 are charged leptons with different flavours and h is the scalar particle with mass 125. 25 GeV, in a two-Higgs-doublet model where all the Yukawa-coupling matrices conserve the lepton flavours but the Majorana mass terms of the right-handed neutrinos break the flavour lepton numbers. We find that (1) the decays ± 1 → ± 2 γ require large Yukawa couplings and very light right-handed neutrinos in order to be visible, (2) the decays Z → + 1 − 2 will be invisible in all the planned experiments, except in a very restricted range of circumstances, but (3) the decays h → + 1 − 2 may be detected in future experiments for rather relaxed sets of input parameters. The well-established phenomenon of neutrino oscillations [1,2] implies that the family lepton numbers are not unbroken symmetries of Nature. Therefore, other processes that violate those symmetries, like the two-body decays ± 1 → ± 2 γ, Z → + 1 − 2 , and h → + 1 − 2 may occur. (Here, h is the recently discovered scalar particle with mass m h = 125. 25 GeV, and 1 and 2 are charged leptons with different flavours.) In the Standard Model (SM) those decays only appear at the one-loop level and they are suppressed by a GIM-like mechanism [3], due to the lightneutrino masses being very small and almost identical when compared to the Fermi scale. As a consequence, in the SM those lepton-flavour-violating (LFV) decays have very small rates and are, in practice, invisible. This fact renders them all the more inviting to explore both experimentally, as windows to New Physics, and theoretically, in extensions of the SM.
In this paper we numerically compute the above-mentioned decays in a simple extension of the SM. That extension is a particular case of the scheme proposed in Ref. [18], which is characterized by the following features: • It is a multi-Higgs-doublet model.
• It has three right-handed neutrinos (RHν), with Majorana masses that enable a type-I seesaw mechanism.
• All the Yukawa-coupling matrices are diagonal in lepton-flavour space, because of the invariance of the dimension-four terms in the Lagrangian under the lepton-flavour symmetries.
• The violation of the family lepton numbers arises only softly, through the dimension-three Majorana mass terms of the RHν.
In Ref. [18] the above-mentioned decays have been computed analytically within that general scheme. In this paper we check that analytical computation, but express the amplitudes through Passarino-Veltman (PV) functions. That allows us to use the resulting formulas in highprecision numerical computations and to establish the impact of the separate amplitudes on the branching ratios (BRs) of the LFV decays. Although our analytical results allow one to study the LFV decays in a model with an arbitrary number of scalar doublets, in this paper we only perform the numerical computation in the context of a simple version of the two-Higgsdoublet model (2HDM). The branching ratios of the LFV decays predicted by seesaw models like ours are usually small due to the strong suppression from the very large RHν Majorana masses [19,20]. A recent paper [21] about the LFV Higgs decays in the framework of a general type-I seesaw model with mass-insertion approximation concludes that the maximal decay rates are far from the current experimental bounds. The inverse seesaw model, a specific realization of low-scale seesaw models, might yield larger decay rates [22][23][24][25].
There is a large number of theoretical papers on the LFV decays, therefore we refer only to some of them, grouping them according to the decays under consideration, since most of the research has been conducted on individual types of decays: • LFV decays of charged leptons were analyzed in the context of the inverse seesaw model [22], of effective field theory [26][27][28][29], of 2HDMs [28,30,31], and of the flipped 3-3-1 model [32]. The current experimental and theoretical situation for these decays is reviewed in Ref. [33].
As mentioned above, the three types of LFV decays have been analyzed mostly separately, but there are also studies that endeavour to combine all three types together [60,61]. Correlations among separate decay rates may exist, and some LFV decays may be constrained by other LFV decays. Some constraints could appear in several models, while other constraints operate only in specific models. For example, Ref. [24] shows that Z → τ ± µ ∓ is constrained by τ ± → µ ± γ in the inverse seesaw model; constraints on the Z decays from the LFV decays of charged leptons also emerge in the 2HDM [41] and in the minimal 3-3-1 model [42]. The authors of Ref. [48] claim that h → τ ± µ ∓ is constrained by τ ± → µ ± γ in their specific 2HDM, but in Ref. [46] no such constraints have been found in the type-III 2HDM.
In this paper we perform the numerical study of all nine LFV two-body decays (τ ± → µ ± γ, τ ± → e ± γ, µ ± → e ± γ, Z/h → τ ± µ ∓ , Z/h → τ ± e ∓ , and Z/h → µ ± e ∓ ) in the context of the 2HDM with seesaw mechanism and flavour-conserving Yukawa couplings. Our purpose is to see under which circumstances the decay rates might be close to their present experimental upper bounds-namely, whether one has to resort to either very large or very small Yukawa couplings, to a very low mass of the charged scalar of the 2HDM, or to very low RHν masses. We want to elucidate which are the very relevant and the less relevant parameters of that model for the LFV decays.
In section 2 we review both the scalar and the leptonic sectors of our model. Section 3 contains our main numerical results. The findings of this paper are summarized in section 4. The Passarino-Veltman functions relevant for the analytic computation are expounded in appendix A. The full one-loop analytical formulas for the LFV decays in terms of PV functions are collected in appendices B, C, and D. Appendix E makes a digression on the invisible Z decay width and appendix F reviews some literature on lower bounds on the charged-Higgs mass.
We assume that no other scalar fields exist, except possibly SU (2) singlets of charge either 0 or ±1. The neutral fields ϕ 0 k have vacuum expectation values (VEVs) v k √ 2 that may be complex. We use the formalism of Ref. [62], that was further developed in Refs. [18] and [63]. The scalar eigenstates of mass are n charged scalars H + a (a = 1, . . . , n) and m real neutral scalars S 0 b (b = 1, . . . , m), with n ≥ n d and m ≥ 2n d . The fields ϕ + k and ϕ 0 k are superpositions of the eigenstates of mass according to The matrix U is n d × n and the matrix V is n d × m. In general, they are not unitary; however, there are matricesŨ that are n × n unitary and m × m real orthogonal, respectively. The matrices T and R account for the possible presence in the model of charged-scalar SU (2) singlets and of scalar gauge invariants, respectively. The unitarity ofŨ and the orthogonality ofṼ imply By definition, H + 1 := G + and S 0 1 := G 0 are the 'would-be Goldstone bosons'. Hence [63], where In Eq. (6), s w and c w are the sine and the cosine, respectively, of the weak mixing angle, and e is the electric charge of the proton. Clearly, because of Eqs. (5) and (6), Thus, V † V 11 = 1. In Eq. (3),Ũ is unitary andṼ is orthogonal. Hence, because of Eq. (7), the first columns of T and R are identically zero. Therefore the orthogonality ofṼ implies that, for b = 1, Thus, is real. So, V † V 1b = −ix b is imaginary for all b = 1.

Some interactions of the scalars
The parameters x b in Eq. (9) are important because they appear in the interaction of the neutral scalars S 0 b with two W gauge bosons [63]: Another important interaction is the one of a W gauge boson with one neutral scalar and one charged scalar. It is given by [18] Also relevant in this paper is the interaction of a neutral scalar with two charged scalars. We parameterize it as where the coefficients obey λ aa b = λ * a ab because of the Hermiticity of L. Equation (12) corresponds, when either a = 1 or a = 1, to an interaction of the charged would-be Goldstone bosons. The coefficients for those interactions may be shown-either by gauge invariance or indeed through an analysis of the scalar potential [18]-to be where m a is the mass of the charged scalar H + a and m b is the mass of the neutral scalar S 0 b .

Leptonic sector 2.2.1 The matrices U and X
We assume the existence of three right-handed neutrinos ν R , where = e, µ, τ . We assume that the flavour lepton numbers are conserved in the Yukawa Lagrangian of the leptons: All the 2n d matrices Γ k and ∆ k are diagonal by assumption. The charged-lepton mass matrix M and the neutrino Dirac mass matrix M D are respectively. The matrices M and M D are diagonal just as the matrices Γ k and ∆ k , respectively. Without loss of generality, we choose the phases of the fields R in such a way that the diagonal matrix elements of M are real and positive, viz. they are the charged-lepton masses; thus, The neutrino mass terms are where C is the charge-conjugation matrix in Dirac space. The flavour-space matrix M R is non-diagonal and symmetric; it is the sole origin of lepton mixing in this model. There are six physical Majorana neutrino fields ν i = Cν T i (i = 1, . . . , 6). The three ν L and the three ν R are superpositions thereof [18]: where P L := (1 − γ 5 )/ 2 and P R := (1 + γ 5 )/ 2 are the projectors of chirality. The matrices U and X are 3 × 6. The matrix is 6 × 6 and unitary, hence The matrix U 6 diagonalizes the full 6 × 6 neutrino mass matrix as In Eq. (21), the m i (i = 1, . . . , 6) are non-negative real; m i is the mass of the neutrino ν i . From Eq. (21), The matrix M D is diagonal. Therefore, It follows from Eqs. (20) and (22) that Equation (21) implies M T D X * = U * m . Therefore, .

The interactions of the leptons
The charged-current Lagrangian is The neutral-current Lagrangian is where When extracting the Feynman rule for the vertex from line (26b), one must multiply by a factor 2 because the ν i are Majorana fields.
The charged scalars interact with the charged leptons and the neutrinos through The coefficients in Eq. (28) are given by The neutral scalars interact with the charged leptons and with the neutrinos through When extracting the Feynman rule for the vertex from line (30b), one must multiply by a factor 2 because the ν i are Majorana fields. The coefficients in Eq. (30) are given by Notice that f bij = f bji . The reader may now appreciate the practical computation of the amplitudes for ± 1 → ± 2 γ (appendix B), Z → + 1 − 2 (appendix C), and S 0 b → + 1 − 2 (appendix D). Those amplitudes are expressed in terms of the Passarino-Veltman functions defined in appendix A.

Restriction to a two-Higgs-doublet model
In the numerical computations in this paper, we work in the context of a two-Higgs-doublet model without any scalar SU (2) singlets. We use, without loss of generality, the 'Higgs basis', wherein only the first scalar doublet has a VEV, and moreover that VEV is real and positive: In this basis, G + = H + 1 is the charged would-be Goldstone boson and H + = H + 2 is a physical charged scalar. Thus, the matrix U defined through Eq. (2) is the 2 × 2 unit matrix. Moreover, G 0 = S 0 1 is the neutral would-be Goldstone boson, and where T is a real orthogonal 3 × 3 matrix. Without loss of generality, we restrict T 11 , T 21 , and T 31 to be non-negative-this corresponds to a choice for the signs of S 0 2 , S 0 3 , and S 0 4 , respectively. Without loss of generality, we choose the phase of the doublet Φ 2 in such a way that T 12 + iT 13 is real and non-negative; thus, T 12 + iT 13 = 1 − T 2 11 . The matrix V defined through Eq. (2) is given by Then, 11 1 ±iT 31 ∓iT 21 iT 21 ∓iT 31 1 ±iT 11 iT 31 ±iT 21 ∓iT 11 1 We are interested only in S 0 2 , viz. in the index b = 2. Through the definition (9), From now on we will only use x 2 and we will not mention T and its matrix elements again. We use the notation m h for the mass of S 0 2 ; since S 0 2 is supposed to be the scalar discovered at the LHC, m h = 125.25 GeV. We use the notation m H + for the mass of the charged scalar H + . The neutral scalars S 0 3 and S 0 4 are unimportant in this paper. We use the following notation [65] for the matrix elements of Γ 1,2 and ∆ 1,2 : while (Γ 1 ) = √ 2m /v. Clearly, according to Eq. (15) with v k = vδ k1 , From Eq. (24), We use both Eqs. (31) and the definition (27) to derive The scalar S 0 2 couples to pairs of gauge bosons according to the Lagrangian [63] It couples to the τ and µ leptons through the Lagrangian-cf. Eq. (30a)- Experimentalists usually write viz. with factors κ W,Z,τ,µ that parameterize the deviations from the SM. Detailed limits on those factors have been derived from experiment, see for instance Refs. [66][67][68][69]. In our fits we enforce the conditions These conditions constitute quite strong constraints on x 2 and on the Yukawa couplings γ τ and γ µ . Conditions (44b) and (44d) are displayed in Fig. 1. In the experimental papers, for any given decay mode a coupling modifier is defined as κ 2 i = Γ i /Γ i SM , therefore in our analysis we allow for either positive or negative Re κ µ and Re κ τ , as illustrated in Fig. 1. We parameterize the vertex of S 0 2 with two charged scalars through Eq. (12). We already know from Eqs. (13) that The value of λ 222 , i.e. of the coupling H − H + S 0 2 , depends on the scalar potential. If we write the quartic part of the scalar potential of the 2HDM in the standard notation [70] then [64] The coupling λ 222 is important for S 0 there is a diagram for that decay wherein S 0 2 attaches to H − H + . However, in practice that diagram gives amplitudes (D5) that are always much smaller than the dominant amplitudes (D4) and (D11). We have found that, for −1 < λ 3 < 7 and |Re λ 7 | < 1.5 [71], the branching ratios BR S 0 are almost completely independent of λ 222 . 3 Thereafter we have kept λ 3 = Re λ 7 = 1 fixed.

Fit to the lepton-mixing data
The lepton mixing matrix U is in the charged-current Lagrangian (25). It is a 3 × 6 matrix. We must connect it to the standard PMNS 3 × 3 unitary matrix. In order to make this connection we use the seesaw approximation [72][73][74][75][76], which is valid when the eigenvalues of M R are very much larger than the (diagonal) matrix elements of M D . The 3 × 3 symmetric matrix is diagonalized by an unitary matrix V as where the n p (p = 1, 2, 3) are real and positive. It follows from Eqs. (48) and (49) that In our fitting program we input the PMNS matrix V , 4 the Yukawa couplings d e,µ,τ , and the n p . We firstly write the matrix M D given by Eq. (38). We then determine M R through Eq. (50). We use M R and M D to construct the 6 × 6 matrix We diagonalize the matrix (51) through the unitary matrix U 6 as in Eq. (21). We thus find both U , viz. the 3 × 6 upper submatrix of U 6 , and the neutrino masses m i (i = 1, . . . , 6). Because the seesaw approximation is very good, one obtains m i ≈ n i for i = 1, 2, 3 and moreover the 3 × 3 left submatrix of U turns out approximately equal to V . Finally, we order the heavy-neutrino masses as m 4 ≤ m 5 ≤ m 6 . Since the inputted n p are many orders of magnitude below the Fermi scale, the matrix elements of M R are much above the Fermi scale unless the Yukawa couplings d are extremely small. Therefore, when we lower the inputted d , we lower the heavy-neutrino masses.
For the n p we use the light-neutrino masses. The cosmological bound [77] is together with the squared-mass differences ∆ solar = n 2 2 −n 2 1 and ∆ atmospheric = |n 2 3 − n 2 1 |, that are taken from phenomenology. The lightest-neutrino mass is kept free; we let it vary in between 10 −5 eV and ∼ 0.03 eV for normal ordering (n 1 < n 3 ), and in between 10 −5 eV and ∼ 0.015 eV for inverted ordering (n 3 < n 1 ); the upper bound on the lightest-neutrino mass is indirectly value of λ 222 becomes quite relevant. However, in that very contrived case the branching ratios of S 0 2 → + 1 − 2 become very close to zero and, therefore, uninteresting to us, since in this paper we are looking for the possibility of largish LFV branching ratios. 4 Recall that in our model there is conservation of the flavour lepton numbers in the Yukawa couplings and therefore the charged-lepton mass matrix is diagonal from the start. provided by the cosmological bound (52). The smallest n p cannot be allowed to be zero becausê n −1 appears in Eq. (50). For the matrix V we use the parameterization [ where ≡ s 13 exp (iδ), c pq = cos θ pq , and s pq = sin θ pq for (pq) = (12), (13), (23). Three different groups [78][79][80] have derived, from the data provided by various neutrinooscillation experiments, values for the mixing angles θ pq , for the phase δ, and for ∆ solar and ∆ atmospheric . The results of the three groups (especially the 1σ bounds) are different, but in Ref. [78] the values of the observables are in between the bounds of Refs. [79] and [80]. In this paper we use the 3σ data from Ref. [78] Table 2: The neutrino-oscillation parameters used in our fits [78].
α 21 and α 31 are kept free in our analysis.
3 Numerical results

Details of the computation
We have generated the complete set of diagrams for each process in Feynman gauge by using the package FeynMaster [81] (that package combines FeynRules [82,83], QGRAF [84], and FeynCalc [85,86]) with a modified version of the FeynRules Standard-Model file to account for the six neutrinos, for lepton flavour mixing, and for the additional Higgs doublet. The amplitudes generated automatically by FeynMaster were expressed through Passarino-Veltman (PV) functions by using the package FeynCalc and specific functions of FeynMaster. All the amplitudes were checked by performing the computations manually. The results of these computations are presented in Appendices B, C, and D. For numerical calculations we made two separate programs, one with Mathematica and another one with Fortran. Because of the very large differences among the mass scale of the light neutrinos, between 10 −5 eV and 0.1 eV, the mass scale of the charged leptons, between 100 keV and 1 GeV, and the mass scale of the heavy neutrinos, between 100 GeV and 10 16 GeV, there are both numerical instabilities and delicate cancellations in the calculations. These numerical problems could be solved with the high-precision numbers that Mathematica allows. However, this strongly slows down the calculations. Fortunately, numerical inaccuracies occur only for very small values (less than 10 −30 ) of the branching ratios (BRs), therefore we were able to use a program written with Fortran to implement the minimization procedure and to find BRs within ranges relevant to experiment. Some parts of the Fortran code (such as the module for matrix diagonalization) have used quadruple precision to avoid inaccuracies, but most of the code has used just double precision so that the computational speed was sufficient for minimization. The final results were checked with the high precision afforded by Mathematica.
The numerical computation of the PV functions was performed by using the Fortran library Collier [87], which is designed for the numerical evaluation of one-loop scalar and tensor integrals. A major advantage of Collier over the LoopTools package [88] is that it avoids numerical instabilities when the neutrino masses are very large, even when one only uses double precision. The integrals were checked with Mathematica's high-precision numbers and Package-X [89] analytic expressions of one-loop integrals.
In the fits of subsection 3.4, in order to find adequate numerical values for the parameters we have constructed a χ 2 function to be minimized: In Eq. (54), • n is the total number of observables to be fitted; this is usually nine, since we fit the BRs of the nine LFV decays in order to find them within the ranges accessible to experiment.
• Θ is the Heaviside step function.
• O v i is the computed value of each observable.
• O b i is the experimental upper bound on the observable, which is given in table 1.
• k is an appropriately small number that short-circuits the minimization algorithm when The χ 2 function (54) works well even when the calculated BRs and the experimental upper bounds differ by many orders of magnitude. We have performed the fits in subsection 3.4 by minimizing χ 2 with respect to the model parameters-the Yukawa couplings d , δ , and γ , and the PMNS-matrix parameters. The mass of the lightest neutrino and the parameter x 2 were randomly generated before the minimization of the χ 2 function, in order to be able to explore the full range of the neutrino masses and the full range of x 2 . In the fits of subsection 3.4 the mass of the charged scalar H + was usually kept fixed, just as the parameters λ 3 and Re λ 7 of the scalar potential (46). The minimization of χ 2 is not an easy task because of the large number of model parameters that, moreover, may differ by several orders of magnitude, and because there is always a large number of local minima. However, we don't try to find absolute minima, i.e. BRs as close as possible to the experimental upper bound; our purpose is rather to search under which circumstances the decay rates may be in experimentally accessible ranges.
The inputted values of the masses of the leptons and bosons were taken from Ref. [66]. We have used s 2 w = 0.22337 and e ≡ √ 4πα, where α = 1/137.036 is the fine-structure constant. The neutrino-oscillation data are in table 2.
We introduce the shorthands BR( ), BR(Z), and BR(h) for the branching ratios of the We also define the lower bound Y min = 10 −6 and the upper bound Y max = √ 4π ≈ 3.5 on the moduli of the Yukawa coupling constants.

Benchmark points
We produce in table 3 two benchmark points (BPs). For those two BPs the neutrino mass ordering is normal, the neutrino squared-mass differences and the lepton mixing angles take their best-fit values in table 2, the mass of the charged scalar is 750 GeV, and the parameters λ 3 and λ 7 of the scalar potential are both equal to 1. The first nine rows of table 3 contain the inputted values of the Yukawa couplings; the next two lines have the computed values of κ µ and κ τ ; in the next four lines one finds the inputted values of the lightest-neutrino mass m 1 , of the Majorana phases α 21 and α 31 , and of the non-alignment parameter 1 − x 2 ; the next three lines have the computed masses of the heavy neutrinos, ordered as m 4 ≤ m 5 ≤ m 6 ; the last nine lines display the computed branching ratios.
In benchmark point 2 (BP-2) only the BR(h) are sufficiently large to be observed in the future, while the BR( ) and BR(Z) are negligibly small. Benchmark point 1 (BP-1) indicates that very small values of the Yukawa couplings d and large values of the Yukawa couplings δ are required in order to obtain BR( ) in experimentally reachable ranges. BP-2 shows that, if only the BR(h) are accessible, then the Yukawa couplings may all be in the range [0. 1,1]; in that case, since the d are not very small, the heavy-neutrino masses are quite large.

Evolution of BRs
In this subsection we discuss the behaviour of the BRs when we vary some input parameters of the benchmark point 1 of the previous section, while the other input parameters of that point remain fixed.
In order to visualize the impact of the Yukawa coupling constants on the BRs, we have fixed their ratios in the same way as in BP-1, viz. 5 3.6 3.9 Table 3: Two benchmark points. In the third column, the symbol '-' stands for a tiny number 10 −20 . The values of the input parameters absent from the first column are given at the beginning of subsection 3.2.
We change either only d τ , or only δ τ , or only γ τ , and we let the other Yukawa couplings vary together with them through the fixed ratios (55). All the other input parameters keep the values of BP-1.
In Fig. 2 we display the BRs against d τ , while the three δ and the three γ are kept equal to their respective values of BP-1. It should be noted, in the upper and lower horizontal scales of  10 -20  10 -20 a small factor in the decay width, cf. Eq. (C2), the predicted BR(Z) are smaller by more than six orders of magnitude than the present experimental upper bounds.
We observe a completely different behaviour of the BR(h) in the right panels of Fig. 2 (the top panel is a zoom of part of the bottom one): the BR(h) achieve values comparable to the experimental upper bounds for a wide range of d τ , viz. even when the heavy-neutrino masses are quite large.
The main message of Fig. 2 is that all nine BRs would be very small if there were no charged scalars H ± . The contributions to the amplitudes from the diagrams with H ± increase some BRs in some circumstances by several orders of magnitude.
The BRs behave differently when plotted against the Yukawa couplings δ , as shown in Fig. 3. We see that all the BRs increase with increasing absolute value of δ ; the BR( ) and BR(h) become visible in planned experiments when |δ | 2 (for appropriate values of the   In Figs. 2 and 4 we have seen that the behaviour of BR(h) is different from the one of BR(Z) and BR( ). This happens because of different amplitudes, but also because of additional parameters, viz. x 2 and the triple-scalar couplings λ 3 and λ 7 , that arise in the diagram of Fig. 17 where h attaches to two charged scalars with couplings given by Eqs. (45) and (47). However, due to the small factor 1 − x 2 2 in the second term of Eq. (47), the impact of λ 7 on BR(h) is almost imperceptible. On the other hand, λ 3 may have a strong influence on BR(h). This happens only for extremely small values of 1 − x 2 , though; and, for such extremely small values of 1 − x 2 , BR(h) is anyway much too small to be measurable. This is displayed in Fig. 5. In the cases that we are interested in, viz. when the BR(h) are rather large, the exact value of λ 3 is unimportant. For the sake of simplicity, from now one we assume λ 3 = λ 7 = 1 everywhere.

-15
10 -11 10 -7 10 -3 10 -15 10 -12 10 -9 10 -6 10 -3 Excluded by CMS 2018 As shown in Fig. 5, 1 − x 2 has a strong impact on BR(h). It is also important for making BR(h) and BR( ) simultaneously close to the experimental bounds. Indeed, the BR(h) may be made sufficiently large, for a wide range of the Yukawa couplings d and for sufficiently large values of the δ , just by varying 1 − x 2 . The strong impact of 1 − x 2 and of the γ on BR(h) allows one to adjust BR(h), together with BR( ), to be close to the experimental upper bounds-but for a quite restricted range of d and δ , because large BR( ) require extremely small d and rather large δ . If, on the other hand, one attempts to fit only BR(h), then both 1 − x 2 and the Yukawa couplings may be much more relaxed, as shown in BP-2 of table 3.
In Fig. 6 we illustrate the evolution of the BRs when the mass of the charged scalar m H + is changed, while the other parameters are kept fixed at their values of BP-1. One observes that when m H + increases the BRs mostly decrease monotonically as (for the other parameters fixed in their values of BP-1) Eventually, when m H + ∼ 10 7 GeV for BR( ) and BR(Z), and when m H + ∼ 10 9 GeV for BR(h), the BRs settle at their SM values. This illustrates the decoupling of H + . One also sees in Fig. 6 that, for ( 1 , 2 ) = (µ, e), there is near m H + = 750 GeV a partial cancellation of amplitudes that leads to a sudden drop of BR (µ ± → e ± γ); our benchmark point 1 has profited from that effect for attaining BR(µ ± → e ± γ) smaller than its experimental upper bound. In Fig. 7 we display the BRs as functions of the Majorana phase α 31 , with the other input parameters kept fixed at their values of BP-1. Here too, for ( 1 , 2 ) = (µ, e) there is a sudden drop of the branching ratios when α 31 = 1.06, which is precisely the value that we have utilized in benchmark point 1. A similar behaviour of the green lines also occurs with other parameters, besides m H + (Fig. 6) and α 31 (Fig. 7). Hence, the values of the parameters must be chosen very carefully if we want to find all six BR( ) and BR(h) simultaneously close to their experimental upper bounds. The main difficulty arises because the upper bound on BR (µ ± → e ± γ) differs from the upper bound on BR (τ ± → µ ± γ) by five orders of magnitude. Fortunately, our minimization procedure allows this to be done quite efficiently.

Fitting the BRs
In this model there is a large number of input parameters. We have performed a minimization procedure in order to find adequate values for all of them. For each set of input parameters, we have computed the branching ratios of the nine LFV decays; we have then selected sets of input parameters for which all six BR( ) and BR(h) are simultaneously between the current experimental upper bounds and the future experimental sensitivities. 6 Since in this subsection we use a fitting procedure, we must enforce definite bounds on the input parameters, lest they acquire either much too small or much too large values. We adopt the following conditions: • The neutrino-oscillation parameters, viz. the mixing angles θ 12 , θ 13 , and θ 23 , the Dirac phase δ, and the neutrino squared-mass differences, are varied within their respective 3σ ranges taken from Ref. [78] and reproduced in table 2. The Majorana phases α 21 and α 31 are kept free, i.e. we let them vary from 0 to 2π.
• The lightest-neutrino mass is varied in between 10 −5 eV and either ∼ 0.03 eV for normal ordering of the neutrino masses or ∼ 0.015 eV for inverted ordering. The precise upper bound on the lightest neutrino mass is fixed, for each pair of values of ∆ solar and ∆ atmospheric , by the Planck 2018 cosmological upper bound (52).
• The Yukawa coupling constants d , δ , and γ are assumed to be real (either positive or negative). 7 • The moduli of the Yukawa coupling constants are varied between Y min = 10 −6 (which is the order of magnitude of the Yukawa coupling of the electron) and a perturbativity bound Y max = √ 4π ≈ 3.5.
There are experimental and phenomenological constraints on the mass of the charged scalar m H + , as discussed in Appendix F. The numerical study in the previous subsection (see Fig. 6) shows that, when m H + increases, most BRs decrease. Since we attempt to obtain largish BRs, the fitting procedure always tends to produce the lowest m H + in the allowed range. In our fits we have fixed m H + = 750 GeV, in agreement with the lower bounds of recent global fits [90,91]. We have also kept the triple-scalar couplings fixed, viz. λ 3 = λ 7 = 1, because they do not change much the BRs. Finally, we have checked that all the final points meet the 3σ conditions on the Z invisible decay width in Eq. (E6).
In the figures of this subsection we display three different fits: 1. In the first fit (displayed through blue points and called 'NO' from now on), we have assumed normal ordering of the light-neutrino masses.
2. In the second fit (displayed through red points and named 'IO') there is inverted ordering of the light-neutrino masses.
3. In the numerical analysis 8 we have found that most points have |d | close to the lower bound Y min = 10 −6 . Therefore, we have produced a third fit (displayed through green points and labelled 'Y min = 10 −7 ') that has normal ordering like the first one, but where the lower bound on the moduli of the Yukawa couplings is 10 −7 instead of 10 −6 . 9   Fig. 8. Most points in Fig. 8 have negative κ τ . This allows larger BR (h → τ ± µ ∓ ) and BR (h → τ ± e ∓ ). If one only allows positive κ τ , then in NO it is not possible to reach the future sensitivity for BR (h → τ ± µ ∓ ), except if one allows complex Yukawa couplings. On the other hand, in both the IO and Y min = 10 −7 cases it is still possible to get all three BR(h) above their future sensitivities with positive κ µ and κ τ .
We have found that free Majorana phases permit larger BRs for the decays with τ ± . Thus, it is advantageous to fit the Majorana phases instead of fixing them at any pre-assigned values.
In Fig. 9 we display correlation plots of BR( ) and BR(Z). One sees that there is a correla- Figure 9: Correlation plots between BR( ) and BR(Z) for the three fits of Fig. 8; the points and the notation are the same as in that figure.
In Fig. 10 we display histograms of the moduli of the Yukawa couplings for our three fits. In the first row of panels one sees that, in order to get BR( ) in experimentally reachable ranges, our fits always have very small |d | 10 −5 . If we had set Y min much larger than 10 −6 , then it might not have been possible to obtain BR( ) visible in the next generation of experiments. In the fit with relaxed Y min = 10 −7 the distribution of the |d | is more uniform. The second row of Fig. 10 shows that, in all three fits, |δ e,µ,τ | have values close to the allowed upper bound Y max . In the third row one sees that the |γ e,µ,τ | vary in rather wide ranges, from ∼ 10 −4 to Y max . This happens because the parameter 1 − x 2 has a strong impact on BR(h); for smaller values of the γ , larger values of 1 − x 2 still allow BR(h) to reach experimentally reachable ranges.
In Fig. 11 we display the heavy-neutrino masses m 4 , m 5 , and m 6 for our three fits. Because the |d | are always so small in the fits, the heavy-neutrino masses are very small too. Thus, in NO and IO m 4 lies in between ∼ 0.5 TeV and ∼ 2.5 TeV, and in Y min = 10 −7 it may be as small as 45 GeV. 10 The mass m 5 is in between ∼ 1 TeV and ∼ 10 TeV for all cases, and the mass m 6 is in between ∼ 2.5 TeV and ∼ 100 TeV for NO, and ∼ 10 4 TeV for both IO and Y min = 10 −7 . In all three fits, it is found that the mixing angles θ 12,23,13 , the Dirac phase δ, and the Majorana phases α 21 and α 31 may have values anywhere in their ranges.

Single decays
In the previous subsection we have discussed the results that are obtained when all six BR( ) and BR(h) are simultaneously between the current experimental upper bounds and the future experimental sensitivities. Here we consider the case where only one of those six BRs is above the future sensitivity.
We have found that requiring just one BR( ) to be above the future sensitivity still restricts the Yukawa couplings |d | 10 −3 and |δ | 0.5. Then, because of the small |d |, the heavyneutrino mass m 4 is always of order 1 TeV (except in µ ± → e ± γ for which m 4 may be of order 50 TeV).
Requiring just one BR(h) to be above the future sensitivity restricts |δ | 0.1. The |d | and the heavy-neutrino masses do not need to be very small, as one sees for instance in BP-2 of table 3.
In our model the decay Z → µ ± e ∓ might be observed at the FCC-ee collider in a very restricted range of circumstances, viz. with large |δ e,µ | 4, small |d e | 10 −6 and |d µ,τ | 5 × 10 −6 , and small m H + 500 GeV. Moreover, a very precise finetuning is required, wherein  Figure 11: Histograms of the heavy-neutrino masses for the three fits.
the Yukawa couplings are such that on the hand the decay µ ± → e ± γ has a cancellation of amplitudes leading to its BR being below the experimental upper bound, and on the other hand BR (Z → µ ± e ∓ ) still remains a little above the FCC-ee sensitivity. The other two LFV Z decays Z → τ ± e ∓ and Z → τ ± µ ∓ are in our model always much too suppressed to be visible.

Amplitudes
Numerically, we have found that only a few amplitudes have a substantial impact on the BRs.
For ± 1 → ± 2 γ the amplitudes (B14) are dominant. Specifically, a l,H in Eq. (B14a) gives the main impact. Therefore, the approximate decay width is This yields the following approximate formulas for the BRs: BR τ ± → e ± γ ≈ 5.733 × 10 4 |a l,H | 2 , BR µ ± → e ± γ ≈ 2.580 × 10 10 |a l,H | 2 . (58c) The amplitudeā l,H in Eq. (C9a) is the most important one for the BR(Z). 11 Therefore, Hence, The amplitudes for the Higgs decays differ from those for the other decays. The amplitudes from the self-energy-like diagrams of Fig. 16, with the charged scalar H ± , give the strongest 11 Due to the similarities between a l,H in Eq. (B14a) andā l,H in Eq. (C9a), there are correlations between BR ± 1 → ± 2 γ and BR Z → ± 1 ∓ 2 , as already seen in Fig. 9.
impact on the branching ratios. Specifically, the amplitude d rb, 16(a,b) we have In order to check the correctness of the approximate BRs of Eqs. (58), (60), and (62) we have calculated the asymmetry between the exact BRs and the approximate ones, Using the points of case 'NO', these asymmetries are displayed in Fig. 12. One sees that    BR asymmetry 0.1, which means that the approximate formulas are quite accurate. These approximate expressions for the BRs may be very useful for intermediate calculations of the fitting procedure, where the calculations need to be repeated many times, before the final result is calculated by using the exact expressions. This computational trick has saved us a lot of time.

Summary and conclusions
Here we summarize our main findings: • The amplitudes with the charged scalar are crucial in order to obtain BR ± 1 → ± 2 γ and BR h → ± 1 ∓ 2 in experimentally accessible ranges.
• Because the experimental upper bound on BR (µ ± → e ± γ) is very small, it is often necessary to finetune the values of the input parameters of the model so that the largest amplitudes for that specific decay partially cancel among themselves, while the other five decays ± 1 → ± 2 γ and h → ± 1 ∓ 2 remain experimentally visible in the future.
• The decays ± 1 → ± 2 γ necessitate large values of the Yukawa couplings |δ | 1 and extremely small values of the Yukawa couplings |d | 10 −5 in order to be visible. The latter imply a very low seesaw scale, i.e. rather light right-handed neutrinos.
• In our model the decays Z → ± 1 ∓ 2 correlate with the decays ± 1 → ± 2 γ, i.e. they behave similarly as functions of the parameters. Because of this correlation, the experimental upper bounds on BR ± 1 → ± 2 γ imply that the decays Z → ± 1 ∓ 2 will remain invisible in all the planned experiments.
They might be visible in future experiments without the need to choose either very large or very small Yukawa couplings.
• The Majorana phases have a significant impact on the branching ratios of all the decays.
One should refrain from fixing them at any pre-assigned values.
• Both normal and inverted ordering of the light-neutrino masses may yield decay rates of adequate orders of magnitude.
• When the mass of the charged scalar increases, most BRs decrease. Still, for m H + 1.5 TeV the decays ± 1 → ± 2 γ and h → ±

A Passarino-Veltman functions
The relevant Passarino-Veltman (PV) functions are defined in the following way. Let the dimension of space-time be d = 4 − with → 0. We define Then, and Some PV functions in Eqs. (A2) and (A3) have 1/ divergences that are independent of the arguments of the PV functions. Thus, C 00 p 2 , (p − q) 2 , q 2 , A, B, C = 1 2 + finite terms.
All other PV functions in Eqs. (A2) and (A3) converge when → 0. We next introduce specific notations for some PV functions that are used in appendices B, C, and D. Thus, where m 1 and m 2 are the masses of the charged leptons ± 1 and ± 2 , respectively, m i and m j are the masses of the neutrinos ν i and ν j , respectively, m a and m a are the masses of the charged scalars H ± a and H ± a , respectively, and m W is the mass of the gauge bosons W ± .
We compute the process − 1 (p 1 ) → − 2 (p 2 ) γ (q), where q = p 1 − p 2 . Obviously, If the outgoing photon is physical, then q 2 = 0; but we keep q 2 = 0 for generality. The amplitude for a photon with polarization σ is where S has been defined in Eq. (A1) and e is the electric charge of the proton. Clearly, If T σ in Eq. (B2) is multiplied by q σ and then Eqs. (B1) and (B3) are utilized, one must obtain zero because of gauge invariance. Thus, We have used Eqs. (B4)-that hold even when q 2 = 0-as a check on our calculations. The decay width is, in the rest frame of the decaying − 1 , 12 In our model each of the coefficients a l , . . . , c r is the sum of two contributions, viz.
a l = a l,H + a l,W , . . . , c r = c r,H + c r,W .
The contributions with sub-index H arise from the diagrams in Fig. 13 and are given in Eqs. (B14) below, and the contributions with sub-index W come from the diagrams in Fig. 14 and are given in Eqs. (B22) below. Notice that Fig. 13 includes diagrams with the charged Goldstone bosons G ± ≡ H ± 1 . In all our calculations we utilize Feynman's gauge. Let m a denote the mass of H ± a ; for a = 1 one must use m a=1 = m W because we are in Feynman's gauge.
The charged scalars H ± a couple to the charged leptons and the neutrinos according to Eq. (28). The charged scalars include as a particular case the charged Goldstone bosons. For G ± = H ± 1 , one has [18] where U is the lepton mixing matrix and s w is the sine of the weak mixing angle. Figure 13: The three diagrams for 1 The diagrams in Figs. 13(a) and 13(b) produce where Notice that in our model 12 Instead of Eq. (B6) there is another way to express the decay width, viz.
This agrees with Eq. (7) of Ref. [93], that has a factor e 2 missing, though.

B.2 W ±
Besides the diagrams exclusively with G ± , there are five diagrams with W ± , cf. Fig. 14 Figure 14: The five diagrams for 1 − → 2 − γ with a loop containing W ± .
ures 14(d) and 14(e) have the outgoing photon attaching to W ± G ∓ . Those diagrams produce where the f functions have been defined in Eqs. (A6). The diagram of Fig. 14(c) produces A crucial property of the lepton mixing matrix U in our model is cf. Eq. (20). In spite of f 00 containing a divergence 1/(2 ), a l,c in Eq. (B18a) is finite because of Eq. (B19). Finally, there are the diagrams of Figs. 14(a) and 14(b), producing a l,ab = e 2 s 2 a r,ab = e 2 m 2 m 1 s 2 b l,ab = b r,ab = c l,ab = c r,ab = 0, Thus, the sum total of the diagrams of Fig. 14 is We compute the process Z (q) → + 1 (p 1 ) − 2 (−p 2 ), where q 2 = m 2 Z and Eqs. (B1) hold. The amplitude for a Z with polarization σ is written The decay width in the rest frame of the decaying Z is where and We define so that the coupling of the Z to the charged leptons is given by cf. Eq. (26a). Notice that We shall write the coefficientsā l , . . . ,c r as the sum of three pieces, viz.
We recover the diagrams of Fig. 13, with the photon substituted by a Z. Diagrams 13(a) and 13(b) produce the result in Eq. (B9) with the transformations P L → t l P L and P R → t r P R . Diagram 13(c) produces the result in Eq. (B13) multiplied by t l . Thus, the full result of Fig. 13 with Z instead of γ is [y ia (e 0 + 2e 1 ) + m 2 z ia (e 2 + 2e 12 ) + m 1 w ia (e 1 + 2e 11 )] , We consider the diagrams of Fig. 14 with a Z instead of the γ. They producē a l,W =ā l,de +ā l,c +ā l,ab , . . . (C10) c r,W =c r,de +c r,c +c r,ab .

C.3 Diagrams with two neutrino internal lines
There are also diagrams where the Z boson attaches to the neutrino line as depicted in Fig. 15.
The relevant vertex is given in Eq. (26b). We have +m 2 y ija [m j q ij g 2 + m i q ji (g 0 + g 2 )] +z ija q ij −2g 00 + q 2 g 12 − m 2 2 (g 2 + g 12 + g 22 ) − m 2 1 (g 1 + g 12 + g 11 ) −z ija q ji m i m j g 0 − w ija q ji m 2 m 1 (g 0 + g 1 + g 2 )} , {−m 2 x ija [m i q ij (g 0 + g 2 ) + m j q ji g 2 ] + m i m j w ija q ij g 0 −m 1 y ija [m j q ij (g 0 + g 1 ) + m i q ji g 1 ] + m 2 m 1 z ija q ij (g 0 + g 1 + g 2 ) +w ija q ji 2g 00 − q 2 g 12 + m 2 2 (g 2 + g 12 + g 22 ) + m 2 1 (g 1 + g 12 + g 11 ) , (C15b) b l,15(a) = 1 c w s w n a=1 6 i,j=1 [m i y ija q ji g 1 + m 2 z ija q ij g 12 − m 1 w ija q ji (g 1 + g 11 In Eqs. (C15), and the g functions are defined in Eqs. (A7). Note that the divergences cancel out inā l,15(a) andā r,15(a) . Indeed, because the unitarity of the matrix U 6 of Eq. (19) implies U X T = 0 3×3 ; and 6 i,j=1 We compute the process S 0 b (q) → + 1 (p 1 ) − 2 (−p 2 ), where S 0 b is a physical neutral scalar, i.e. b = 1. Equations (B1) hold and q 2 = m 2 b . The amplitude is written where S was defined in Eq. (A1). The decay width in the rest frame of S 0 b is where D.1 Diagrams in which S 0 b attaches to charged leptons There are self-energy-like diagrams with a loop of either H ± a -diagrams (a) and (b) in Fig. 16or W ± -diagrams (c) and (d) in Fig. 16. The vertex of S 0 b with the charged leptons is given Figure 16: The four self-energy-like diagrams for by Eq. (30a). One obtains where j 0,1,2 are defined by Eqs. (A9). Note that the diagram of Fig. 17 implicitly contains the cases where either H ± a or H ± a (or both) coincide with the charged Goldstone bosons G ± := H ± 1 . In those cases one must use m a=1 = m W together with Eqs. (13).

D.3 Diagrams in which S 0 b attaches to W bosons
We next compute the diagrams in Fig. 18. The relevant terms of the Lagrangian are the ones −L a 2 i 4 k 00 + m 2 1 (k 11 + k 12 + k 1 ) + m 2 2 (k 22 + k 12 + 2 k 1 + 2 k 2 ) −q 2 (k 12 + 2 k 1 ) , −L * a 1 i 4 l 00 + m 2 1 (l 11 + l 12 + 2 l 1 + 2 l 2 ) + m 2 2 (l 22 + l 12 + l 2 ) −q 2 (l 12 + 2 l 2 ) , with the l functions defined in Eqs. (A11). Equations (D6) and (D7) contain no divergences because vanishes if = , since the matrices Γ k are diagonal. Equations (D6) and (D7) include the particular case where a = 1; then, H ± a coincides with the Goldstone bosons G ± . In that particular case one must use m a=1 = m W together with Eqs. (B8) and where x b is the real number defined in Eq. (9). In order to compute diagram 18(c) one must know the vertex of a neutral scalar with two W ± gauge bosons, which is given by Eq. (10) [63]. One then obtains where The neutral scalar S 0 b may also attach to two internal neutrino lines. The relevant diagrams are displayed in Fig. 19. The vertex of the neutral scalars with the neutrinos is given by Eq. (30b). {−x ija f bji m 2 m 1 (g 0 + g 1 + g 2 ) − y ija f bji m i m j g 0 −y ija f * bji 4g 00 + m 2 1 (g 1 + g 11 + g 12 ) + m 2 2 (g 2 + g 22 + g 12 ) − q 2 g 12 +z ija f bji m i m 2 (g 0 + g 2 ) + z ija f * bji m j m 2 g 2 +w ija f * bji m i m 1 g 1 + w ija f bji m j m 1 (g 0 + g 1 ) , −y ija f * bji m 2 m 1 (g 0 + g 1 + g 2 ) − x ija f * bji m i m j g 0 −x ija f bji 4g 00 + m 2 1 (g 1 + g 11 + g 12 ) + m 2 2 (g 2 + g 22 + g 12 ) − q 2 g 12 +z ija f * bji m j m 1 (g 0 + g 1 ) + z ija f bji m i m 1 g 1 +w ija f bji m j m 2 g 2 + w ija f * bji m i m 2 (g 0 + g 2 ) , effectively constrain the branching ratios of LFV processes in our model. This is distinct from LFV studies in the inverse seesaw model [24] or in the effective field theory of the seesaw [39]. That happens because in our case the masses m 4,5,6 of the heavy neutrinos are sufficiently high that the Z can never decay into a heavy neutrino plus a light neutrino, except for very small values |d | 10 −7 of the Yukawa couplings; and because the non-unitarity of the matrix U has a very weak impact on the couplings of the active neutrinos to the Z boson in our model.

F Constraints on the mass of the charged scalar
Direct constraints on m H + may be obtained from collider experiments on the production and decay of on-shell charged Higgs bosons. The search sensitivity is limited by the kinematic reach of experiments, but collider constraints have the advantage of being robust and modelindependent. The bound obtained from direct searches at LEP for any value of tan β is m H + > 78.6 GeV at 95% CL [97]. Combining data of the four LEP experiments, a limit of m H + 80 GeV is obtained [66,98], while m H + 150 GeV may be derived from the searches at LHC [99,100]. Stronger mass limits on m H + may be obtained for specific regions of tan β. Some constraints on m H + from flavour physics depend strongly on the 2HDM Yukawa type, while others are type-independent. Among the flavour processes, the constraints from b → sγ are most stringent due to the constructive interference of the H ± contribution with the SM contribution. For a type-II 2HDM, the lower limit m H + > 480 GeV at 95% CL [101] includes NNLO QCD corrections and is rather independent of tan β. In a recent study [102], the branching ratio of b → sγ enforces m H + 580 GeV at 95% CL both for the type-II and for the flipped 2HDM.
The recent global fits in Refs. [90,91,100,103] give bounds on the charged-Higgs mass for various 2HDM Yukawa types. In those studies only 2HDMs with a Z 2 -symmetric potential are considered, but one may suppose that the bounds would be similar for the general 2HDM. In Ref. [100] it is found that, for the type-II 2HDM, flavour-physics observables impose a lower bound m H + 600 GeV that is independent of tan β when tan β > 1 but increases to m H + 650 GeV when tan β < 1. In Ref. [90], m H + > 740 GeV in both the type-II and flipped 2HDMs, but m H + 460 GeV for the lepton-specific 2HDM. In Ref. [91] on finds m H + 500 GeV or m H + 750 GeV in the aligned 2HDM, depending on the fitted mass range. However, for the type-I and lepton-specific 2HDMs the restrictions on m H + from flavour constraints are weaker [100,103,104].