Neutrinoless double beta decay and the muonium-to-antimuonium transition in models with a doubly charged scalar

The lepton number and flavor violations are important possible ingredients of the lepton physics. The neutrinoless double beta decay and the transition of the muonium into antimuonuim are related to those violations. The former can give us an essential part of fundamental physics, and there are plenty of experimental attempts to observe the process. The latter has also been one of the attractive phenomena, and the experimental bound will be updated in planned experiments at new high-intensity muon beamlines. In models with a doubly charged scalar, not only can those two processes be induced, but also the active neutrino masses can be induced radiatively. The flavor violating decays of the charged leptons constrain the flavor parameters of the models. We study how the muonium-to-antimuonium transition rate can be as large as the current experimental bound, and we insist that the updated bound of the transition rate will be useful to distinguish the models to generate the neutrinoless double beta decay via the doubly charged scalar. We also study a possible extension of the model to the left-right model which can induce the neutrinoless double beta decay.


Introduction
The neutrinoless double beta decay (0ν2β decay) is one of the key phenomena to probe lepton number violation. The experimental searches for the 0ν2β decay from 76 Ge [1], 130 Te [2], and 136 Xe [3] are ongoing, and the 0ν2β half-value periods of them are bounded to be more than about 10 25 -10 26 years depending on the nuclei. The next-generation experiments including other isotopes are planned [4,5,6,7,8,9,10], and hundred times higher sensitivity of the half-lives is targeted. Together with precise data from the neutrino oscillation experiments and cosmological observations of the sum of the neutrino masses, it is expected to advance our understanding of the mass structure of neutrinos.
On the lepton physics, the high intensity muon beamline is being upgraded [11,12], and precise measurements of the muon properties are planned. The transition of the muonium (µ + e − ) into antimuonium (µ − e + ) (Mu-to-Mu transition) [13,14,15] is one of the phenomena to search at the facilities. It has been pointed that the transition and the 0ν2β decay can have a close relationship [16]. The Mu-to-Mu transition violates lepton flavor numbers (∆L µ = −∆L e = 2). Although lepton flavor violation (LFV) can be easily considered in models beyond the standard model (SM), the LFV decays, such as µ → eγ [17] and µ → 3e [18], and µ-e conversion in nuclei [19] give severe constraints to the models. Those LFV processes violate lepton flavor numbers with single units (∆L e , ∆L µ = ±1). The Mu-to-Mu transition can be allowed theoretically even if the LFV decays of the charged leptons are severely suppressed. The Mu-to-Mu transition experiment [20] provides the bound of the coefficients of the fourfermion operators to induce the transition (∼ 10 −3 in the unit of the Fermi coupling constant G F ), and the bound will be lowered to O(10 −4 )G F by the experiments at the facilities with the upgraded muon beamlines [21,22,23].
The (1,1) element of the neutrino mass matrix, m ee , gives a contribution of the 0ν2β amplitude from the ordinary diagram with the long-range exchange of the Majorana neutrino. If there are particles beyond the SM and there is a new lepton number violating interaction, the 0ν2β amplitude can be generated [24,25] even if m ee is small. The doubly charged scalar is one of the representative new particles that can induce the 0ν2β decay via short-range interactions. In the models with the doubly charged particle, the active neutrino mass can be also induced radiatively (see Ref. [26] for a review of the radiative neutrino mass models). The Majorana masses of the active neutrinos and the 0ν2β decay are linked by the lepton number violation. Various possibilities of such types of the models that can induce the 0ν2β decay have been considered [27,28,29,30,31,32,33,34]. The doubly charged scalar couplings to the leptons can carry lepton flavor numbers, and therefore, the Mu-to-Mu transition can be also induced via the doubly charged scalar exchange [35,36] in those models. The Mu-to-Mu transition can be a tool to select models to induce the 0ν2β decay, which is what we focus on in this paper. If SU(2) L singlet scalars, h + and k ++ , with their hypercharges Y = 1 and Y = 2 are added in the two-Higgs-doublet model (2HDM), the 0ν2β decay can be induced via the diagram given in Fig.1. We call that diagram a cocktail diagram [37]. In the model, the charged scalar η − in the SU(2) L Higgs doublet is mixed with the singlet h + , and the mixing violates the lepton number. The active neutrino masses can be induced at the two-loop level. We call this model the Liu-Gu model since they discuss both the 0ν2β decay and the loop-induced neutrino mass matrix in their model in Ref. [31]. The half-value period of the 0ν2β decay via the cocktail diagram can be as large as the current experimental bound. The singlet h + can have bilinear couplings to the left-handed lepton doublets ℓ, and then, two-loop induced neutrino masses can be obtained. The model with the ℓℓh + coupling is called the Zee-Babu model [38,39,40,41]. In the original Zee-Babu model, only one Higgs doublet is enough to generate the neutrino masses and thus, the 0ν2β decay is induced by m ee from the ordinary diagram. As a consequence, the half-value period of the decay is much larger than the current bound in the normal neutrino mass hierarchy. Authors have shown that the Mu-to-Mu transition rate can be as large as the current experimental bound in the Zee-Babu model [42]. In the Liu-Gu model, on the other hand, the Mu-to-Mu transition rate is bounded to be much smaller than the current bound due to experimental constraints of the LFV decays of the charged leptons. Therefore, the 0ν2β decay and the Mu-to-Mu transition are useful to distinguish those two models.
Can both the 0ν2β decay amplitude and the Mu-to-Mu transition rate be as large as the current bounds? In this paper, we show that the Mu-to-Mu transition rate can be as large as the current bound avoiding the constraints from the LFV decays and the 0ν2β decay is induced by the cocktail diagram in Fig.1 if the ℓℓh + coupling is added in the Liu-Gu model (in other words, the η + -h + mixing is simply added in the Zee-Babu model). Those amplitudes are proportional to the e R e R k ++ coupling, and the scalar spectrum is predictable to induce the half-value period of the 0ν2β decay and the Mu-to-Mu transition. In the Liu-Gu model, the type II (or type X) 2HDM is preferable to induce the neutrino masses with suppressing the LFV decays. In order to induce the sizable Mu-to-Mu transition operator, it is possible to apply the type I (or type Y) 2HDM in our hybrid model, and the scalars in the doublets can be lighter while satisfying the experimental constraints from their decays. If the scalars in the doublets are lighter, the half-value period of the 0ν2β decay can be shorter. We will also describe the 0ν2β decay and the Mu-to-Mu transition in the left-right model in which the h + and k ++ scalars can be unified in a SU(2) R triplet. This paper is organized as follows: In Section 2, we introduce the 2HDM with the SU(2) L singlet scalars h + and k ++ , and we present an overview of the Liu-Gu and Zee-Babu models. In Section 3, we describe the 0ν2β decay amplitude via the cocktail diagram in Fig.1, and show the numerical calculation of the half-value periods. In Section 4, the loop-induced neutrino masses are given. In Section 5, the Mu-to-Mu transition in Liu-Gu and Zee-Babu models are discussed. In Section 6, we exhibit the analyses for our model that can make both amplitudes of the 0ν2β decay and Mu-to-Mu transition as large as the current experimental bounds, and discuss how the experimental data can distinguish the models. In Section 7, we develop the model to the left-right model with a SU(2) R triplet. Section 8 gives the conclusion and discussion of this paper. In Appendix A, the scalar spectrum in the 2HDM with the h + and k ++ scalars are given. The condition to avoid the charge symmetry breaking vacua is also noted. In Appendix B, the loop functions for the cocktail and box diagrams are given. In Appendix C, the loop functions for the neutrino masses are shown. In Appendix D, the LFV constraints are summarized. In Appendix E, the scalar spectrum in the left-right model with a SU(2) R triplet is given.

Overview
We introduce two SU(2) L doublets, Φ 1 and Φ 2 , and two SU(2) L singlets, h + and k ++ , with hypercharges Y = 1 and Y = 2, respectively. The doublets acquire vacuum expectation values (vevs) to break SU(2) L × U(1) Y : and with a notation: c β = cos β and s β = sin β, so that the vev of the doublet η is zero. The components in Φ and η are and ω − and ω 0 are the would-be Nambu-Goldstone bosons eaten by the W and Z gauge bosons.
To avoid flavor changing neutral currents, a discrete symmetry is often assumed in the 2HDM [43,44]. For the Yukawa coupling of the up-type quarks, the one with Φ 2 is conventionally chosen. For the Yukawa couplings of the down-type quarks and the charged leptons, four types of choices are possible and they are referred to as type I, II, X, and Y. Leaving out the selection for the down-type quarks, Table 1 shows the choices of the discrete symmetry for the charged leptons. The coupling of the charged scalar in η to the leptons are where M e is the charged lepton mass matrix, M e = diag(m e , m µ , m τ ). The lepton couplings to h + and k ++ are given as follows: where ℓ denotes the left-handed lepton doublets, ℓ = (ν L , e L ) T , and ℓ c i ℓ j (with a proper SU(2) L contraction) can be written as The coupling f ij is anti-symmetric under the flavor indices, while the coupling g ij is symmetric. The scalar potential in terms of Φ 1 , Φ 2 , h + , and k ++ is (2.10)  where only those allowed by the Z 2 symmetry are written for the quartic couplings in terms of Φ 1 and Φ 2 . We define the vev-dependent masses of h + and k ++ as The neutral scalars Φ 0 and η 0 are mixed (with the mixing angle β − α), and their mass eigenstates are denoted by h 0 and H 0 . The masses of the neutral scalars, h 0 , H 0 , and A 0 , are denoted as m h , m H , and m A . The SM(-like) Higgs boson is h 0 , and m h = 125 GeV. The singly charged scalars, η + and h + , can be mixed by the Φ 1 Φ 2 h + = Φηh + term, and we denote the mass eigenstates as h + 1 and h + 2 , and their masses as m h + 1 and m h + 2 . We denote the mixing angle of them as φ, and The expressions of the scalar masses are given in Appendix A. We choose the Z 2 -charge of k ++ to be + in order to allow the g coupling. For the Z 2 -charge assignment of h + , we consider two choices as given in Table 1. They correspond to the following models: 1. Z 2 -charge A: Liu-Gu model [31].
The ℓℓh + coupling is forbidden. The Φ 1 Φ 2 h + = Φηh + coupling is allowed. The cocktail diagram in Fig.1 induces the 0ν2β decay whose amplitude can be sizable to observe the process in the near future. Neutrino masses can be generated by the diagram in Fig.2. We note that h − k ++ Φ 1 Φ 2 coupling is allowed and the 0ν2β amplitude via the cocktail diagram and the neutrino masses can be also generated via the h − k ++ Φ 1 Φ 2 coupling in addition to the h + h + k −− coupling. The loop-induced neutrino mass matrix is (2.13) The Mu-to-Mu transition rate is restricted due to the LFV constraints, as we will explain later.
The ℓℓh + coupling is allowed. The Φ 1 Φ 2 h + coupling is forbidden 1 . The cocktail diagram in Fig.1 is absent and the 0ν2β decay is induced only by the usual contribution from m ee . Neutrino masses can be generated by the diagram in Fig.3, and (2.14) The Mu-to-Mu transition rate can be sizable to observe in the near-future experiments.
It is important that the 0ν2β decay and Mu-to-Mu transition are available to select the models with the Z-charges A and B since only one of their amplitudes can be sizable. The Φ 1 Φ 2 h + scalar trilinear coupling has a dimension, and thus, the trilinear coupling can be a soft-breaking term of Z 2 -charge B. If the soft-breaking term is added in the model of Z 2 -charge B, it is possible that both amplitudes of the 0ν2β decay and Mu-to-Mu transition are large. The purpose of this paper is to investigate the possibility that both can be large enough to be observed in the near future.

Neutrinoless double beta decay via the cocktail diagram
In this section, we describe the 0ν2β contribution from the cocktail diagram [29] in Fig.1, and we show the numerical calculation.
In order to execute the triangle loop calculation, we remark that the part of the propagator for the neutral scalars in the mass eigenstates can be written as (3.1) If cos(β − α) = 0 and m H = m A , the contribution to the 0ν2β decay amplitude vanishes. Since we will work in the alignment limit cos(β − α) → 0 to avoid experimental constraints, we ignore the h 0 /H 0 part and only consider the H 0 /A 0 part in the following.
The operator induced by Fig.1 is (3.2) and the coefficient is where the loop function F t is given in Appendix B, and The current bounds of the 0ν2β half-value periods by GERDA [1], CUORE [2], and Kamland-Zen [3] are We use numerical values of nuclear matrix elements (NME) given in Ref. [45]. The half-value period for 136 Xe gives a little stronger bound to the amplitude than those of 76 Ge and 130 Te due to the NME, and we obtain the numerical bound in the case of λ hkΦΦ = 0 as follows: The experimental constraints for the doubly charged scalar mass are found in [46,47]. In the following, we suppose m k ∼ 1 TeV, which can be allowed by the experiment.
In Fig.4, we plot the half-value period for Xe as a function of m + h2 for m 2 k /m hhk = 1 TeV, g ee = 0.3, and various masses of the scalars, m H , m A , and m h + We suppose λ hkΦΦ = 0. As described in Appendix A, the scalar trilinear coupling should be roughly smaller than the quadratic masses to avoid charge symmetry breaking, and we assume the η + -h + mixing to plot the figure as We note that the amplitude vanishes when m h + 1 = m h + 2 . It is expected that the 0ν2β decay can be observed in the future if |g ee | ∼ O(0.1) and the masses of H, A, and h + 1 (those come from the doublets, roughly) are less than several hundreds of GeV.
We note that a 0ν2β operator, (u L d R ) 2 (e R e c R ), can be generated at the tree level via a u L d R η + coupling in the Liu-Gu model. However, the u L d R η + coupling is related to the down quark Yukawa coupling, and the tree-level contribution to the 0ν2β amplitude is less than the loop contribution. In general 2HDM (so-called the type III model), the u L d R η + coupling can be sizable, and the tree-level contribution can be made large compared to the 2HDM with the Z 2 charge, though we do not touch it in this paper.

Loop-induced neutrino mass matrices
The neutrino masses are induced by the diagrams given in Fig.2 and Fig.3 in the Liu-Gu and Zee-Babu models, respectively, and the flavor part of the induced neutrino mass matrices are given in Eqs. (2.13) and (2.14). In this section, we show more explicit expressions of the neutrino mass matrices and give the solutions in terms of the f and g coupling matrices to reproduce the observed neutrino oscillations.
The η − couplings to the leptons are given as In the Liu-Gu model, the loop-induced neutrino mass matrix by the diagram in Fig.2 is In the Zee-Babu model without h + -η + mixing, the neutrino mass matrix induced by the diagram in Fig.3 is The loop functions L 1 and L 2 are given in Appendix C. It can be seen that the flavor-dependent parts are given as in Eqs.(2.13) and (2.14), and the model parameters determine the proportional coefficients. We note that a three-loop diagram with a cocktail diagram [37] can also induce the neutrino mass in the Liu-Gu model. The flavor structure is also given in Eq.(2.13).
In the Liu-Gu model, the solution of the g coupling is simply In the Zee-Babu model, we need a little explanation to show its solution for g and f matrices. First of all, the f coupling is an antisymmetric matrix, and thus, the rank of the neutrino mass matrix is two, and the mass matrix (in the normal mass hierarchy) is given as where U is the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) neutrino mixing unitary matrix, and m 2 , m 3 are neutrino masses (with a Majorana phase in a convention, arg(m 2 /m 3 ) = α 2 ). We write the mixing matrix as where u i 's are column vectors. We set From the orthogonality of the column vectors, one obtains (4.10) We find that gives the algebraic solution of f M e g * ZB M e f T ∝ M ν [42]. Due to f u 1 = 0, a i 's are free. The flavor structure of the g coupling in the Zee-Babu model can be controlled by the three free parameters. The solution in the inverted mass hierarchy (m 3 = 0) is obtained similarly by choosing f ij = f 0 ǫ ijk (u 3 ) k . We consider only the case of the normal mass hierarchy in this paper.
It is clear from this algebraic solution that all the neutrino oscillation parameters, three mixings, a Dirac phase δ, two masses (m 1 = 0), and one Majorana phase (α 2 ), can be realized in the Zee-Babu model in principle. The (1,1) element of the neutrino mass matrix, which is an important parameter for the usual light neutrino exchange diagram for the 0ν2β decay, is thus simply calculated in the Zee-Babu model as and we obtain |m ee | ≃ 1-4 meV, (4.13) by using the current experimental measurement of the oscillation parameters [48].

Mu-to-Mu transition and LFV constraints to realize the neutrino matrices
Either the f coupling or the off-diagonal elements of the g coupling matrix can induce the LFV decay processes. In this section, we review the LFV constraints in the Liu-Gu and Zee-Babu models. Our concern is the Mu-to-Mu transition from the exchange of the doubly charged scalar k ++ . We describe whether the transition rate can be as large as the current experimental bound, avoiding the LFV constraints.
There are five independent four-fermion Mu-to-Mu transition operators [49]. The induced operator by the k ++ exchange is and The bound from the transition experiment [20] is The LFV constraints are a bottleneck to obtaining a sizable Mu-to-Mu transition rate if the g coupling matrix has off-diagonal elements. Appendix D summarizes the LFV constraints. To satisfy the stringent constraint from the µ → 3e decay, g eµ has to be tiny. The constraint from µ → eγ decay gives a constraint to the size of not only g eµ but also g eτ (if g µτ = 0). In the case of the Liu-Gu model in Eq.(4.6), the (1,3) element of the neutrino mass matrix M ν has to be non-zero to reproduce the 12 and 13 neutrino mixings since g eµ is tiny. Then, g eτ needs to be sizable, and consequently, |G 2 |/G F is bounded to be less than O(10 −5 ) due to τ → 3e in the Liu-Gu model (see Eq.(D.13) and Ref. [42]).
In the Zee-Babu model, on the other hand, both g eµ and g eτ can be made to be zero using the free parameters a i in Eq.(4.11). As a result, |G 2 |/G F ∼ O(10 −3 ) is possible. The description of the LFV constraints in the Zee-Babu model is found in the literature [50,51,52,53,54]. The algebraic solution of the Zee-Babu model is given in Eqs.(4.9) and (4.11). Since the elements of the f coupling are determined by the first column of the PMNS unitarity matrix, the factor f 0 of the f coupling matrix needs to be small due to the µ → eγ constraint given in Eq.(D.3). Then, the elements of the g coupling matrix need to be sizable to reproduce the proper neutrino masses for a scalar mass spectrum m h + 1 , m k ∼ 1 TeV. Due to m e ≪ m µ ≪ m τ , the g ee element is easily more than unity. To avoid g ee ≫ 1, we need to use one of the degrees of freedom of the parameters a i . Since we need to make g eµ → 0, one of g eτ and g µτ cannot be adjusted. The restriction from the τ + → e − µ + µ + decay to g eτ element is more stringent due to the charged lepton mass hierarchy. Consequently, g µτ cannot be made small, and the µ → eγ decay from the f coupling and the τ → 3µ decay from the g coupling restrict the mass parameters in the Zee-Babu model. The Mu-to-Mu transition rate is bounded by the τ → 3µ and τ + → µ − e + e + decays, and the bound of the transition rate is just below the current experimental bound [42].
6 Analyses for our model As we have explained, the cocktail diagram can generate the 0ν2β amplitude and the 0ν2β half-value periods of Ge and Xe can be just above the current bound in the Liu-Gu model. In the Zee-Babu model, on the other hand, the 0ν2β amplitude from the cocktail diagram is not generated and the half-value period is 1000-10000 times larger than the current bound since Φηh + coupling is absent and the 0ν2β decay is induced only by m ee given in Eq.(4.13). The Mu-to-Mu transition rate is far below the current bound in the Liu-Gu model due to the LFV constraints, while it can be just below the current bound in the Zee-Babu model. These are useful to distinguish the models. Even in the model with Z 2 -charge B in Table 1 for the Zee-Babu model, the Φ 1 Φ 2 h + term can be employed as a soft-breaking term in principle. Then, the cocktail diagram can generate the 0ν2β amplitude, and the half-value period can be just above the current bound if g ee is sizable. At that time, the Mu-to-Mu transition rate can be as large as the current bound satisfying the LFV constraints, which we will verify in this section. 2 The two-loop induced neutrino mass matrix in the Z 2 -charge B model with the soft-breaking where L LG and L ZB are given in Appendix C. The dimensionless coupling λ hkΦΦ is forbidden by the Z 2 charge. This equation is rewritten as We comment that there can be a three-loop contribution in the Liu-Gu model, and the parameter R can be redefined by containing the contribution in Eq.(6.2).
We assume g eµ = g eτ = 0 to satisfy the LFV bound and fix the value of g ee to give a contribution to the 0ν2β decay by the cocktail diagram. Then, there are three parameters in g and three parameters in f coupling matrices. Equation (6.2) contains six equations. Therefore, for fixed values of R, L ZB , m 2 k /m hhk , and neutrino oscillation parameters in M ν , one can solve the equations for g µµ , g µτ , g τ τ , f eµ , f eτ , and f µτ . The solutions are screened by the LFV constraints, especially for |f µτ f eτ | from the µ → eγ bound. From the assumption that g eµ = g eτ = 0, the solar neutrino mixing and 13-mixing have to be generated from f . Therefore, roughly saying, R should be smaller than O(0.1) to satisfy the µ → eγ bound.
To consider the LFV constraints of the model parameters, the size of g µµ is important. Roughly speaking, the numerical size of g µµ is determined by the (2,2) element of Eq.(6.2) as To reproduce the atmospheric neutrino mixing, one needs As a result, the τ + → µ − e + e + decay process constrains the solution of a smaller R for fixed m k . The size of G 2 /G F is constrained to be less than O(10 −3 ) from the τ + → µ − e + e + and τ → 3µ decay processes. See Eqs.(D.14) and (D.16). In Fig.5, we show a plot of G 2 /G F as a function of R to depict the statements above. The solid line shows the plot in the case of the Liu-Gu model (f = 0 in Eq.(6.2)). The (1,2) element of M ν is chosen to be zero (namely g eµ = 0 to satisfy the stringent µ → 3e constraint), which means that the neutrino oscillation parameters are related (see Ref. [42]). The ratio of g eτ /g µµ is fixed as a solution to reproduce M ν . We choose g ee = 0.3 and m k = 1 TeV for , it violates the bound of τ → 3e process, and therefore, we need to enlarge m k to satisfy the bound. Then, the solid line becomes flat for R < ∼ O(1), which can be understood by Eq.(D.13) for fixed g eτ /g µµ . We note that the µ → eγ decay can restrict G 2 /G F by Eq.(D.5) for smaller values of R. The circles are plotted as the solutions of Eq.(6.2) that satisfies the LFV constraints for g ee = 0.3, m h + 2 = 1 TeV, L ZB = 2/(16π 2 ) 2 , and m hhk = m k = 1 TeV. In order to obtain the points, the neutrino mass and mixing parameters [48], as well as phases, are randomly generated within experimental errors. The asterisks are plotted as the solutions similarly for g ee = 0.3, m h +  In those solutions, the branching ratios of the µ → eγ decay are smaller than the experimental bound in our setup of the singly charged scalar mass spectrum due to |f eτ f * µτ | < 0.001 (see Appendix D).
In our model of Z 2 -charge B with Φ 1 Φ 2 h + coupling, both the Mu-to-Mu transition rate and the 0ν2β amplitudes can be large for R ≃ O(0.001)-O(0.1). This size of R can be realized if r < 1, which corresponds to the type I/Y 2HDM. In fact, the experimental constraints of the scalar decays to fermions can allow light scalars in the doublets [55], especially in the type I 2HDM. On the other hand, the Liu-Gu model with Z 2 -charge A needs a larger r to reproduce the neutrino masses while suppressing the LFV decays for m k ∼ m hhk ∼ 1 TeV, and then, the light scalars are constrained from the experimental data of the scalar decays to fermions. The lighter scalars can produce the smaller half-value periods as shown in Fig.4.

Extension to the left-right model
In the left-right model whose gauge symmetry is SU(3) c × SU(2) L × SU(2) R × U(1) B−L , the singly and doubly charged scalars can be unified in the SU(2) R triplet ∆ (with B − L = 2), The two SU(2) L doublets are unified into a bidoublet, (2, 2, 0) under SU(2) L × SU (2) [71,72,73]. As a result, the amplitude from the one-loop cocktail diagram can be larger than the amplitudes at the tree-level W R exchanges. In this section, we discuss those contributions with a new box loop contribution in the model of the extension to the left-right model. Before going to the left-right model, we simply add a singlet fermion N to the 2HDM with h + and k ++ . The coupling and mass of N are given as The κ coupling can induce the 0ν2β operator from the box diagram in Fig.6. The coefficient of the operator in Eq.(3.2) is and the loop function F b is given in Appendix B. We comment that one can also consider the Dirac neutrino coupling, ℓNΦ. However, the coupling can generate too large neutrino masses (even if it only couples to the inert doublet, ℓNη, and the Dirac neutrino mass is absent [74]), and therefore, the coupling should be too small to contribute to the 0ν2β amplitude. In Fig.7, we plot the ratio of the amplitudes from the cocktail and box diagrams. We choose m 2 k /m hhk = 1 TeV and m h + 2 = 1 TeV in the plot. The box loop can give a subleading contribution to the 0ν2β amplitude. If m hhk /m 2 k ≪ 1/(1 TeV), the cocktail contribution is suppressed and the box contribution can become the leading contribution to the decay amplitude. Now we consider the left-right model. The right-handed charged lepton and the SM singlet fermion N are contained in a SU(2) R doublet ℓ R , and the g coupling in Eq.(2.8) is extended as As we have mentioned at the beginning of this section, there are tree-level contributions from the W R gauge boson to the following 0ν2β operator: Let us estimate the tree-level contributions for the current experimental bound of M WR > ∼ 5 TeV. Since our purpose here is to compare the tree-level contributions with the one-loop contributions, we here describe the contributions only from the diagrams in Figs.8 and 9. Surely, there are additional contributions that pick W L -W R mixing. The contributions from the diagrams are We note that the W − R -W − R -∆ ++ term is given by the gauge interaction where W R = τ a 2 W a R and τ a 's are the Pauli matrices. We also note that m N = 2g ee ∆ 0 if additional singlet fermions are not introduced and the heavy neutrino mixings are absent. We obtain where A 0 R is the coefficient of the operator that gives the current experimental bound of the 0ν2β half-value period of Xe by using the NME for the operator given in Ref. [45]. We note that M WR ≃ g R ∆ 0 if the vev of ∆ 0 dominantly breaks SU(2) R × U(1) B−L . In the minimal model where the vev of ∆ 0 is the only source to break SU(2) R × U(1) B−L → U(1) Y , ∆ + is absorbed to the W R boson. The η − -∆ − mixing φ is determined only by the vevs irrespective of the detail of the scalar potential. We obtain See Appendix E for more details. As a consequence, the contribution from the cocktail diagram via W L boson exchanges cannot overcome the contribution from the tree-level diagram in Fig.9. If one introduces two SU(2) R triplets that have vevs to break SU(2) R × U(1) B−L or an additional SU(2) R doublet (1, 2, 1/2) under (SU(2) L , SU(2) R , (B − L)/2), the η − -∆ − mixing depends on the detail of the scalar couplings. For example, in the inverse seesaw model [75] which is often considered to obtain the tiny neutrino masses more naturally in the TeV-scale left-right model, the additional SU(2) R doublet is motivated to be introduced. If the vev of the SU(2) R doublet mainly generates the W R boson mass and the vev of ∆ 0 is smaller than it, the η − -∆ − mixing can become a free parameter (independent from the gauge boson mass ratio) in the non-minimal model. Then, as in the 2HDM with h + and k ++ , the amplitude of the one-loop contribution from the cocktail diagram can be as large as the current bound, and can overcome the contribution from the tree-level W R exchanges. The box loop (Fig.6) (with κ = √ 2g ee ) can also overcome the contribution from the right-handed neutrino in Fig.8 if m N is heavier than several hundred GeV. We note that there is a W + L -W − R -η 0 term from the gauge interaction and the W R boson can propagate in the loops in the cocktail and box diagrams, though those  contributions to the 0ν2β decay amplitude are small due to the heaviness of the W R boson. Surely, there are also tree-level diagrams of the heavy Majorana neutrino N exchanges with active-sterile neutrino mixings [59,60,61,62,63,64,65,66,67,68,69]. However, the activesterile mixings need to be enlarged to obtain the sizable amplitudes. The large active-sterile mixing can induce the active neutrino masses at the loop level [76], and then, the Majorana mass of N should not be large to obtain the sizable 0ν2β amplitude in those cases to utilize the large active-sterile mixings [77] and the box loop contribution can be more significant for In the left-right model, the tiny active neutrino masses can be generated at the tree level even if the g coupling matrix is diagonal and the LFV constraints are avoided. Therefore, the Mu-to-Mu transition operator via the tree-level ∆ ++ exchange can be as large as the current experimental bound if the doubly charged scalar mass is about 1 TeV, as discussed in Ref. [42]. The Yukawa interactions of the quarks and leptons, q L : (3, 2, 1, 14) The up-and down-type quarks' mass matrices are linear combinations of Y q1 and Y q2 , and the charged leptons and the Dirac neutrino mass matrices are linear combinations of Y ℓ1 and Y ℓ2 .
In the inverse seesaw model, we introduce three generations of gauge singlet fermions S i and a SU(2) R doublet, H R : (1, 1, 2, 1/2). The interactions, (7.15) give the active neutrino mass matrix as for small elements of µ S , where m D is a Dirac neutrino mass matrix, and m S = d H 0 R . As described, for ∆ 0 ≪ H 0 R , the H + R boson mainly becomes the longitudinal component of W R boson, and the ∆ + and ∆ ++ components in the SU(2) R triplet ∆ respectively behave as h + and k ++ to induce the 0ν2β decay via the cocktail diagram in Fig.1.

Conclusion
The 0ν2β decay and the Mu-to-Mu transition are investigated in the 2HDM with SU(2) L singlet scalars, h + and k ++ . The current bounds of the 0ν2β decay translate to m ee < ∼ 0.1 eV if the ordinary diagram of the Majorana neutrino exchange dominates the decay amplitude. It is well known that the observation of the decay just above the current bounds of the half-lives does not necessarily mean m ee ∼ 0.1 eV in the models beyond the SM. The e R e R k ++ and k ++ h − h − interactions, as well as the mixing between h + and a charged scalar in the SU(2) L Higgs doublet, can violate the lepton number. The 0ν2β decay can be induced via the diagram given in Fig.1, and the half-life can be as large as the current bound even if m ee is small. We have shown that both the 0ν2β decay and the Mu-to-Mu transition rate can be as large as the current experimental bounds in the general setup of the 2HDM with h + and k ++ that induces the neutrino mass matrix by two-loop diagrams in Figs.2 and 3, while the experimental constraints from the LFV decays of the charged leptons are satisfied. The search for both the 0ν2β decay and the Mu-to-Mu transition will be important to select models and specify the parameter space of the model.
If we consider another possibility that the Mu-to-Mu transition is induced in connection with neutrinos, the type-II seesaw model immediately comes to mind. The doubly charged scalar is contained in the SU(2) L triplet ∆ L for the type-II seesaw term, and ℓℓ∆ L terms can induce an operator for the Mu-to-Mu transition. However, the constraints of the LFV decays require neutrino mass degeneracy (even in the case of the normal mass hierarchy) to obtain the sizable Mu-to-Mu transition operator. Consequently, the cosmological bound of a total of the neutrino masses [78] gives a limitation on the Mu-to-Mu transition rate (see Refs. [22,42]). If the type-I seesaw contribution is added, both the 0ν2β decay amplitude and the Mu-to-Mu transition rate can be as large as the current bound. However, one needs a cancellation in m ee between the type-I and type-II seesaw contributions, which impairs the predictivity of the model. In the model we propose in this paper, the amplitudes of the 0ν2β decay and the Mu-to-Mu transition are proportional to the e R e R k ++ coupling (without a cancellation), and the scalar spectrum is more predictable to induce the half-lives of the 0ν2β decay and the Mu-to-Mu transition. The experimental results at the LHC and ILC will give us information on the possible scalar spectrum [79,80,81,82,83,84,85,86]. This paper claims that the Mu-to-Mu transition, which will be updated in the near-future experiments, can be one of the players to sort out those models which induce the 0ν2β decay.

Note added
We correct the discussion in Section 6. The induced neutrino mass M ν should include not only the two-loop contributions described by Eq. (6.1), but also a one-loop one. Even after the correction, our main claim is unchanged that the measurement of the Mu-to-Mu transition is important to sort out the models which induce the 0ν2β decay. As discussed in Section 2, the f ij ℓ i ℓ j h + coupling is forbidden by the Z 2 -charge A (given in Table 1), and only "Liu-Gu" contribution in Eq. (2.13) is generated by the two-loop diagram in Fig. 2. On the other hand, the Φ 1 Φ 2 h + scalar trilinear term is forbidden by the Z 2 -charge B, and only "Zee-Babu" contribution in Eq. (2.14) is generated by the two-loop diagram in Fig. 3. The loop-induced effective neutrino mass operators are (ℓΦ † 1 )(ℓΦ † 1 ) or (ℓΦ † 2 )(ℓΦ † 2 ). We have introduced a breaking term of the Z 2 discrete symmetry and mix the Liu-Gu and Zee-Babu contributions comparably to see if both the Mu-to-Mu transitions and 0ν2β decay width can be large enough to be observed by near-future experiments. However, when the breaking term is introduced, an effective (ℓΦ † 1 )(ℓΦ † 2 ) operator can be generated by an one-loop diagram shown in Fig. 10, as discussed in the original Zee model, If the two-loop Liu-Gu and Zee-Babu contributions are comparable, the one-loop Zee contribution dominates the neutrino mass matrix. As it is known, the original Zee model alone cannot reproduce the neutrino oscillation parameters. The proper treatment for our purpose is, therefore, that we introduce a small f eτ coupling to the Liu-Gu model; the two-loop Liu-Gu and the one-loop Zee contributions are comparable (the two-loop Zee-Babu contribution is tiny). Then, the (1,3) element of M ν is generated by the one-loop Zee contribution, and the observed neutrino oscillation parameters can be reproduced even in the case of g eµ = g eτ = 0 to avoid too large LFV decay rates. Therefore, the Mu-to-Mu transition can become large, as discussed in the text. The situation is much simpler to solve the equation given in Eq. (6.2).
Thus, the third and subsequent paragraphs in Section 6 should be replaced with the following sentences: The induced neutrino mass in the Z 2 -charge B model with the soft-breaking where L LG is given in Appendix C. The dimensionless coupling λ hkΦΦ is forbidden by the Z 2 charge, and M ZB ν is ignored because the second order of f is negligible unless s 2φ is extremely small.
We assume g eµ = g eτ = 0 to satisfy the LFV bound and fix the value of g ee to give a contribution to the 0ν2β decay by the cocktail diagram. We can realize such a suitable size of g ee by adjusting the lightest neutrino mass without conflict with the constraints for neutrino masses. Since Eq. (8.2) is a linear equation for any elements of g and f , we can simply solve for the couplings. We find two trivial solutions: 1. We assume f µτ = 0. Then, we obtain for g µµ , g τ τ , and g µτ , and for f eµ and f eτ . When the neutrino mass matrix is given, the couplings are numerically determined.
In Fig. 11, we show a plot of |G 2 | /G F as a function of r 2 s 2 2φ . The red solid line shows the plot in the case of the Liu-Gu model (f = 0 in Eq. (8.2)). The (1,2) element of M ν is chosen to be zero (namely g eµ = 0 to satisfy the stringent µ → 3e constraint), which means that the neutrino oscillation parameters are related (see Ref. [42]). The ratio of g eτ /g µµ is fixed as a solution to reproduce M ν . We choose g ee = 0.3 and m k = 1 TeV for r 2 s 2 2φ > ∼ O (1). For r 2 s 2 2φ < ∼ O(1), it violates the bound of τ → 3e process, and therefore, we need to enlarge m k to satisfy the bound. Then, the red solid line becomes flat for r 2 s 2 2φ < ∼ O(1), which can be understood by Eq. (D.13) for fixed g eτ /g µµ . We note that the µ → eγ decay can restrict G 2 /G F by Eq. (D.5) for smaller values of r 2 s 2 2φ . The blue solid line as the solution of Eq. (8.2) where g ee = 0.3, L LG = 0.1/ 16π 2 2 , and m hhk = 1 TeV. We set m k = 1 TeV for r 2 s 2 2φ > ∼ O(0.1), while we choose m k (> 1 TeV) to satisfy the τ + → µ − e + e + and τ → 3µ bounds for r 2 s 2 2φ > ∼ O(0.1).
Since |f eτ f µτ | is sufficiently small, the constraint from µ → eγ does not affect the plot. In order to obtain the line, we use the best-fit values of the neutrino mass and mixing parameters [48]. In those solutions on the blue solid line, the 0ν2β contribution from the cocktail diagram can be a little above the current experimental bound if the scalars in the doublets are light as shown in Section 3. We find that the model parameters for r 2 s 2 2φ ≃ 0.1 can induce the Mu-to-Mu transition with a detectable transition rate just below the current bound shown by the horizontal-dotted line.
2. We assume g µτ = 0. We obtain the same expressions as Eqs. (8.3) and (8.4) for g µµ , g τ τ , f eµ , f eτ , and f µτ . Since the assumption of g µτ = 0 suppresses the LFV τ decays, we can more easily evade the constraints from LFV searches.
The mass eigenstates and the mixing angle are conventionally defined as and the so-called alignment limit (cos(β − α) → 0) without decoupling is obtained if The CP-odd neutral scalar mass is The singly charged scalar mass term is where When the scalar trilinear coupling, such as m hhk and m ΦΦh in Eq.(2.10), is large, charge symmetry breaking vacua can become a true vacuum as often discussed in supersymmetric models [87]. The trilinear coupling, therefore, needs to be bounded to avoid the phase transition to the charge symmetry breaking vacua. For example, for a direction of h + = k ++ , the value of the potential V is positive at the charge breaking vacua and the charge breaking is not a true vacuum if This is a sufficient condition since it is acceptable that the electroweak symmetry breaking vacua is metastable in the cosmological time scale. As a rule of thumb, the trilinear coupling should be smaller than the order of the quadratic masses.

B Loop functions for neutrinoless double beta decay
What we need to calculate loop functions is , we obtain the triangle loop part of the cocktail diagram as follows: , .

(C.2)
We define loop functions for the respective diagrams as In Fig.12, we plot the loop functions for D LFV constraints for f and g coupling matrices The couplings f and g in Eq.(2.8) can induce LFV decay processes. Let us summarize the constraints from the processes. The branching ratio of the µ → eγ decay induced by the f coupling is The experimental bound of the branching ratio is [17] Br(µ → eγ) < 4.2 × 10 −13 , and the f eτ f * µτ product is constrained as The branching ratio of the decay induced by the g coupling is Br(µ → eγ) = 3α 16π 4(g ee g * µe + g eµ g * µµ + g eτ g * µτ ) 3G F m 2 k 2 , (D. 4) and the experimental bound of the coupling can be similarly written. Since the g coupling can also induce the Mu-to-Mu transition operator in Eq.(5.1), it is convenient to express the bound as follows: |G 2 | G F < 3.8 × 10 −6 g ee g µµ g ee g * µe + g eµ g * µµ + g eτ g * µτ . (D.5) The g coupling can induce the µ-e conversion in nuclei [19] with a large ln(m 2 µ /m 2 k ) correction [88]. The current experiment gives a weaker constraint than the µ → eγ bound.
The τ → µ(e)γ decays [89], though the bounds are not stringent. We comment that the h + -η + mixing does not induce a new one-loop bound from the µ → eγ decay, since an active neutrino mass is inserted in the loop diagram with the h + -η + mixing. In the two-loop level, there can be a diagram given in Fig.13. We can verify that the two-loop contribution does not provide a significant constraint as long as the f couplings satisfy the one-loop level constraint and the g µτ coupling is suppressed rather than g µµ to reproduce the neutrino mass matrix.