A model of neutrino mass and dark matter with large neutrinoless double beta decay

We propose a model where neutrino masses are generated at three loop order but neutrinoless double beta decay occurs at one loop. Thus we can have large neutrinoless double beta decay observable in the future experiments even when the neutrino masses are very small. The model receives strong constraints from the neutrino data and lepton flavor violating decays, which substantially reduces the number of free parameters. Our model also opens up the possibility of having several new scalars below the TeV regime, which can be explored at the collider experiments. Additionally, our model also has an unbroken $Z_2$ symmetry which allows us to identify a viable Dark Matter candidate.


Introduction
Almost all the extensions of the Standard Model (SM) directed towards an explanation for the neutrino masses brings in the possibility of lepton number violation (LNV) as an outcome. It is well known that neutrinoless double beta decay (0νββ decay) which is a convincing signature for LNV, will be an inevitable consequence if the neutrino has Majorana mass. If the main contribution to 0νββ decay proceeds through the Majorana neutrino propagator, depending on the spectrum of the neutrino masses, the expected rate for 0νββ decay might be too small to be observed in the experiments. But there exist scenarios where the dominant mechanism for 0νββ decay is not controlled by the Majorana neutrino propagator. In such cases we can have the possibility of large 0νββ decay even when the neutrino Majorana masses are small. Many studies have been performed in this direction in the past (see Refs. [1][2][3] for a general overview, Refs.  for specific models 4 and Refs. [28][29][30][31] for effective field theory (EFT) approaches). In Ref. [31], the authors performed an EFT analysis of the different ways of generating 0νββ decay and light neutrino masses by including operators involving only leptons, Higgs and gauge bosons. This led to a class of interesting models where 0νββ decay was generated at tree level whereas neutrino masses would appear only at two-loops (see Refs. [32] for example models in this category).
The model in Ref. [32] contains an SU (2) L singlet doubly charged scalar like in the Zee-Babu model [33][34][35], an SU (2) L triplet scalar with hypercharge +1 and a real singlet scalar. A Z 2 symmetry, which is later broken spontaneously, is required to prevent tree-level neutrino masses. The model is economical in the sense that it contains no new fermions and by design, it gives new contributions to 0νββ decay, which, in principle, can be large. Additionally, it has a rich phenomenology which can be probed through the searches for the lepton flavor violating (LFV) signals and/or the direct searches for the new scalars in the collider experiments.
In this article we will present a simple variation of the model in Ref. [32]. Our new model will have the same field content as in Ref. [32], except that the Z 2 symmetry will not be broken spontaneously. Consequently, 0νββ decay will now occur at one-loop whereas neutrino masses will appear at three-loop order. The fact the Z 2 is exact makes the model simpler and allows for a viable Dark Matter (DM) candidate: the lightest of the electrically neutral Z 2 -odd particles. On the other hand, the model keeps all the virtues of the previous model: very predictive neutrino mass matrix, large 0νββ decay decay, rich lepton flavour violation phenomenology and new scalars which are in the sub-TeV region and therefore, are within the reach of the collider experiments in the near future.
Our paper will be organized as follows. In Sec. 2 we lay out the scalar field content and the physical spectrum of our model. In Sec. 3 we discuss the 0νββ decay and the bounds that follow from it. Neutrino masses and constraints from LFV decays are discussed in Sec. 4 and Sec. 5 respectively. We analyze the feasibility of DM in Sec. 6. Finally, we summarize our findings in Sec. 7.

The model
The scalar sector of our model contains the following fields: where, the numbers inside the curly brackets associated with the fields represent their transformations properties under SU (2) L and U (1) Y respectively. The normalization for the hypercharge is such that the electric charges of the component fields are given by, Q = T 3 + Y . The fields, χ and σ are odd under an additional Z 2 symmetry which has been introduced to prevent the occurrence of tree-level neutrino masses as well as to ensure the stability of the DM particle. The most general scalar potential involving these fields is given below: where 'Tr' represents the trace over 2 × 2 matrices and Φ = iσ 2 Φ * , with σ 2 being the second Pauli matrix. We can take all the parameters in the potential to be real without any loss of generality.
For the leptonic Yukawa sector, we have the following Lagrangian: where, L L = (ν , ) T L denotes the left-handed lepton doublet and R represents the right-handed charged lepton singlet. C is the charge conjugation operator. We choose to work in the mass basis of the charged leptons which means, Y e is a diagonal matrix with positive entries and f is a complex symmetric matrix with three unphysical phases.

The scalar spectrum
We do not want to break the Z 2 symmetry spontaneously. Denoting by v the vacuum expectation values (vev) of the doublet the minimization conditions read  Table 1: Z 2 parity assignments to the physical particles in our model.
After spontaneous symmetry breaking (SSB) we represent the doublet and the triplet as follows: where, ω and ζ represent the Goldstones associated with the W and Z bosons respectively. Because of the unbroken Z 2 symmetry, only h t and σ can have nontrivial mixing. This leads to a very simple scalar spectrum as described below.
The masses for the doubly charged particles are given by, The mass of the singly charged scalar is given by, The pseudoscalar mass is given by, From Eqs. (6), (7) and (8) it is easy to see that the following correlation holds: In the CP even sector, the SM-like Higgs arises purely from the doublet, Φ, with mass m 2 h = 2λ Φ v 2 . For the other two Z 2 -odd scalars, we obtain the following mass matrix: This mass matrix can be diagonalized by the following orthogonal rotation: with, m 2 and, tan 2α = 2B where we have implicitly assumed that 'S' is the lighter mass eigenstate. One can easily find the following relations: which imply, Combining Eqs. (11) and (13c) we can express λ 6 in terms of the physical parameter as follows: The splittings between different scalar masses can be constrained further from the electroweak T -parameter. The expression for the new physics contribution to the T -parameter is given by where, θ W and M W are the weak mixing angle and the W -boson mass respectively. The function, F (m 2 1 , m 2 2 ), is given by, Taking the new physics contribution to the T -parameter as [36] ∆T = 0.05 ± 0.12 , we will require our model value of the T -parameter to be within the 2σ uncertainty range. For small sin α, this leads to |m H − m χ ++ | 100 GeV.
In passing, combining Eqs. (9) and (14), we note that two types of scalar mass hierarchies are possible depending on the sign of λ Φχ , or, In both cases, m κ ++ can be arbitrary in principle.

Estimation of 0νββ decay
For new scalar masses of O (1 TeV), the Majorana mass matrix element, M ee , will be very small (see Sec. 4 for details). As a result, the usual neutrino exchange diagram will contribute negligibly to 0νββ decay. The main contribution to the 0νββ decay amplitude has been displayed in Fig. 1. From the diagram in Fig. 1 we can easily estimate the effectiveēe c (ūd) 2 interaction giving rise to 0νββ decay where I β is a dimensionless function of the scalar masses running in the loop which is expected to be O (1). For illustration, we have chosen the common scale of the loop to be the mass of the pseudoscalar part from Figure 1: One-loop diagram, in the mass insertion approach, contributing to neutrinoless double beta decay.
the scalar triplet, m A . Of course the diagram in Fig. 1 is only one of the contributions in the mass insertion approach which allows us to give an estimate. A complete calculation of the function I β in the physical basis has been presented in Appendix A yielding values for I β which are slightly smaller than one in the range of masses of interest, I β ∼ 0.1. We will use these values for our estimates.
On the other hand, upcoming experiments are expected to be sensitive to lifetimes of order 10 27 -10 28 yrs [42], i.e. a reduction factor on the coupling of about one order of magnitude. Thus, for 0νββ decay mediated by heavy particles to be observable in the next round of experiments we should have 3 4 × 10 −10 . Therefore in order to escape the current experimental bounds but at the same time to entertain the possibility of observing 0νββ decay in the near future, we require 3 to be within the following range: With f ee , λ 6 ≈ 1, µ κ ≈ m A ≈ m κ ++ ≈ 1 TeV and I β ∼ 0.1 we obtain, from Eq. (22), 3 ∼ 10 −9 which falls naturally within the range given in Eq. (23).

Estimation of the neutrino masses
From Eqs. (2) and (3) it is obvious that simultaneous nonzero values for Y e , f ab , µ κ and λ 6 will prevent us from assigning consistent lepton numbers to all the scalar and lepton fields. Therefore, lepton number is broken explicitly and Majorana neutrino masses will be unavoidable. The sample diagram of Fig. 2, in the mass insertion approach, clearly depicts the involvement of all these couplings in a multiplicative manner. Thus, we can parametrize the neutrino mass matrix as follows: where m a denotes the mass of the charged lepton, a , and I ν represents the loop function expected to be of O (1). Detailed expression of I ν in terms of the scalar masses has been presented in Appendix B. Eq. (24) has a very particular and predictive structure, specific for this class of models, which can be constrasted with the observed spectrum of neutrino masses and mixings (see for instance Refs. [31,32,43]).
As before, taking f τ τ , λ 6 ≈ 1 and µ κ ≈ m κ ++ ≈ 1 TeV and I ν ∼ 1 we obtain the following values for the different elements But of course, some of the f ab s can be much smaller than 1. However, not all of the elements of the f matrix are arbitrary as some of them will be constrained from LFV processes. We will discuss these constraints in Sec. 5. But for now we wish to emphasize that the product |f * ee f eµ | will receive strong bounds from µ → 3e as the latter can proceed at the tree-level mediated by κ ++ . Then, one should naturally expect the following hierarchy among the mass matrix elements: which, obviously, can only accommodate a normal hierarchy among the neutrino masses. In Ref. [32] it has been shown that the above hierarchy with can successfully reproduce the observed masses and mixings in the neutrino sector with a prediction of sin 2 θ 13 > 0.008. Eq. (27) will imply the following hierarchy among the Yukawa elements: We shall also assume f ee f eµ in such a way that f * ee f eµ is still sufficiently small to keep µ → 3e decay under control but at the same time allowing for the possibility of large 0νββ decay.
From Eqs. (22) and (24) we see that the dimensionless factor, is common to both. In terms of γ, the explicit expression for M τ τ in Eq. (27) reads: As we will see in Sec. 5, the ratio f τ τ /m κ ++ is bounded from LFV processes as f τ τ /m κ ++ 1.4 × 10 −4 TeV −1 . Plugging this into Eq. (30) we obtain the following bound for γ: Having an explicit expression for the neutrino masses we can compare the light neutrino exchange contributions to 0νββ decay with the ones discussed in Sec. 3. In fact, from Eqs. (24) and (22) we can express the neutrino mass matrix element M ee , which controls the ν contributions to 0νββ decay, in terms 3 , which parametrizes the new contributions Then, it is clear that for small enough m A the new contributions will dominate over the neutrino contributions. How small? Since the nuclear matrix elements are different in the two cases we cannot make a direct comparison. However, we can use that the experimental limit T 0νββ 1/2 ( 136 Xe) > 1.07 × 10 26 yrs [41] translates into two equivalent bounds on 3 and M ee when 0νββ decay is dominated by the new contributions or by neutrino masses respectively: which already include the appropriate nuclear matrix elements. Using these results and taking I β ∼ 0.1I ν we obtain that the new contributions will dominate for m A 15 TeV. Therefore, scalar masses must be relatively light, and this could make the model testable at the LHC and/or in LFV processes.
Experimental Data (90% CL) Bounds (90% CL) Bounds assuming Eq. (28)  [44,45]. Limits on the Yukawa couplings of the doubly charged singlet scalars have been taken from Ref. [46]. The constraints in the third column are obtained from those in the second column assuming Eq. (28) holds. The bound in the third column corresponding to µ → eγ has an additional assumption, f eµ ≈ 0.

Constraints from LFV processes
Constraints from LFV processes come mainly from decays of the type ∓ ∓ c ∓ d will be more important because these decays can occur at the tree-level through the exchange of the doubly charged scalar singlet, κ ±± . These processes along with the kinds of constraints they imply have been reviewed in Ref. [46] in the context of the Zee-Babu model (see also Refs. [35,47]). The experimental data has not changed much since then. In the first two columns of Table 2 we have summarized the experimental data and the corresponding constraints on the Yukawa couplings. In the third column of Table 2 we recast the constraints of the second column assuming the validity of Eq. (28). This allows us to express the constraints in more specific forms. For example, using m e f eτ ∼ m τ f τ τ and m 2 µ f µµ ∼ m 2 τ f τ τ , the constraint from τ → eµµ leads to a direct bound on f τ τ as follows: It is also worth mentioning that, using Eq. (28), the limit from τ → 3e translates into As mentioned earlier, we want to have f ee relatively large to have appreciable 0νββ decay rate in the future experiments. Then we will need f eµ to be vanishingly small to keep the constraints from µ → 3e under control.
Note that, for f ee ∼ O (1) and sub TeV κ ++ , Eq. (35) will imply a stronger bound on f τ τ than Eq. (34). Figure 3: Regions corresponding to the observed relic abundance [48] in the m S -λ S plane for different values of sin α. We have chosen m H = m χ ++ = m κ ++ = 800 GeV as a benchmark for this plot. Current [49,50] and future [51] bounds from direct detection experiments are also marked appropriately.

Dark Matter
Our model has a Z 2 symmetry which remains unbroken after the SSB. Consequently, the particle spectrum can be divided into Z 2 -even and odd sectors as shown in Table 1. Among the Z 2 odd neutral scalars, S, being the lightest, is a promising candidate for DM. Notice that S is and admixture of the real singlet and the triplet, and therefore, it will feel both, Higgs and gauge interactions 5 . In spite of that, one can parametrize its couplings with the SM-like Higgs boson as follows: with, λ S = 1 2 2λ Φσ cos 2 α − 2 √ 2λ 6 sin α cos α + (λ Φχ + λ Φχ ) sin 2 α .
5 For recent studies of a DM candidate which is an admixture of a scalar singlet and a Y=0 triplet see for instance [52,53].  Table 3: Benchmark values for the input parameters (first row) and other relevant quantities derived from these inputs (second row).
In Fig. 3 we have displayed regions in the m S -λ S plane, which can reproduce the observed DM relic density [48]. For this plot, we have assumed m H = m χ ++ = m κ ++ = 800 GeV and used the MicrOMEGAs package [54] to compute the DM abundance. Note that, the region labeled as sin α = 0 corresponds to the pure Higgs portal scenario. Barring the small window near the Higgs-pole (m S ≈ m h /2, not shown explicitly in the plot), in this case, we need m S 350 GeV [55,56] to evade the direct search bound. It is worth mentioning that in the case of pure Higgs portal, for our choice of benchmark, the DM annihilates through ff , W W , ZZ and hh mainly. All these annihilation channels except hh can only proceed through s-channel h exchange. But as sin α is turned on, we allow for a direct SSV V (V = W, Z) with strength proportional to g 2 sin 2 α. For our choice of positive values for λ S , the new contact diagram will interfere constructively with the h mediated s-channel diagram. 6 This will enhance the annihilation rate for SS → V V once the corresponding threshold is reached. Therefore, we would require lower values of λ S , compared to the pure Higgs portal case, to reproduce the relic abundance. These features have been depicted in Fig. 3 where we can see that a small value of sin α is sufficient to accommodate DM with mass as low as 200 GeV, which can either be discovered or ruled out in the next run of direct detection experiments.

Results and conclusions
Since κ ±± couples directly to the charged leptons, it will be strongly constrained from the same sign dilepton searches at the LHC. Depending on the preferred decay channel of κ ±± , the bound can be as strong as m κ ++ 500 GeV [59,60]. On the other hand, to keep the T -parameter under control, for small sin α, we will need |m H − m χ ++ | 100 GeV (see Eq. (16)). All these considerations together justify our choice of benchmark for Fig. 3. Now, to satisfy Eq. (31) we need to have a large splitting between m H and m S . Keeping these things in mind, we have chosen the first row in Table 3 as a benchmark for the input parameters. Some relevant output quantities that follow from these inputs have also been displayed in the second row of the same table. From the numbers of Table 3 one can easily check that the constraints of Eqs. (23) and (31) and all the bounds in Table 2 are satisfied. Moreover, using Eq. (28) suitable values for f eτ , f µµ and f µτ can be found so that the hierarchy of Eq. (27) is satisfied.
The model has many phenomenological implications that make it special and distinguishable from similar models. To exemplify one such feature, we note that the requirement, M ee , M eµ M eτ , M µµ , M µτ , M τ τ , and consequently NH among the neutrino masses, results in a strong correlation between δ, the CP violating phase of the Pontecorvo-Maki-Nakagawa-Sakata mixing matrix, and the other mixing parameters. For instance, in Fig. 4 we have displayed the allowed region in the plane s 2 23 -δ obtained by the NuFIT collaboration (version 3.2 of 2018) [61,62] (the different coloured contours are 68.27%, 90%, 95.45%, 99% and 99.73% C.L. regions respectively). On top of it we superimpose the correlation obtained from the requirement M ee = M eµ = 0 for the central values of the rest of the mixing parameters (brown dashed line) and the band obtained when they are varied in 1σ. As we can see, the prediction of the model agrees well with the fit, although with some trend to lower values of s 2 23 and δ. Moreover the model also predicts the smallest neutrino mass to be around m 1 ∼ 5 × 10 −3 eV and the two Majorana phases α 1 ∼ 360 • − δ ∼ 130 • and α 2 ∼ α 1 + 180 • ∼ 310 • 7 Figure 4: The NuFIT results [61,62] for the global fit to neutrino data (coloured contours correspond to 68.27% 90% 95.45% 99% 99.73% C.L. regions in the s 2 23 -δ plane) against the prediction of the model for central values of the rest of the mixing parameters (brown dashed line) and the band obtained when they are varied in 1σ.
Eq. (24) allows us to write the couplings f ab in terms of the neutrino masses and mixings up to a global factor. Since these couplings control all the LFV decays mediated by the double charged scalars, all the LFV processes are, in principle, predicted in terms of neutrino masses and mixing parameters which are fixed in our model.
As can be seen from the value of 3 in Table 3, our model opens up the interesting possibility of detecting 0νββ decay in the next generation of experiments even if M ee ∼ 0, but, in addition, is important to remark that the process is quite different from the standard one in which two left-handed electrons are produced. If 0νββ decay is found and proceeds as in the mechanism suggested in this paper, the produced electrons will be right-handed and, therefore, it will be possible, in principle, to distinguish this mechanism by measuring the polarization of the emitted electrons.
We have also found a DM candidate which can reproduce the observed relic abundance yet can survive the current constraints from the direct detection experiments. Furthermore, our model provides the prospect of detecting new scalars with masses below O (TeV) in collider experiments (for LHC studies on lepton number violating singly and doubly charged scalars see for instance [63,64]). Among these new particles, χ ± and χ ±± being Z 2 -odd, cannot decay directly into the SM particles. A search strategy for these kinds of exotic charged scalars can be interesting for the collider studies. Moreover, the decay branching ratios of the singlet doubly charged scalar κ ++ are controlled by the f ab couplings which are fixed in terms of the neutrino mass parameters, therefore, if κ ++ is found at the LHC it will be possible to distinguish this model from other models by comparing the κ ++ leptonic decay branching ratios to neutrino oscillation data and to LFV processes, which also depend on the same couplings.
Appendix A Computation of the loop induced κW W vertex Here we compute the effective κ −− W + µ W µ+ vertex at one loop for vanishing external momenta. Our assumption is justified in view of the fact that the momentum transfers to κ and W -bosons in Fig. 1 are much smaller than the corresponding masses. We write the effective vertex as which, after spontaneous symmetry breaking, emerges from the following gauge invariant operator: After integrating out κ ++ , Eq. (A.2) leads to the following LFV gauge invariant operator [31,32]: We depict in Fig. 5 the three diagrams that contribute to the vertex. Each of these diagrams seem to diverge logaritmically. But one should keep in mind that the neutral scalar exchange must violate lepton number conservation. Thus a large cancellation among the contributions from the three neutral scalars, A, H and S, is expected. After adding all the contributions we obtain an effective neutral scalar propagator of the following form (for Minkowsky momenta) 1 2 where, Φ = v/ √ 2. Evidently, after adding contributions from A, H and S, every diagram in Fig. 5 becomes finite individually. Now we can write the expression of C κW W (defined in Eq. (A.1)) as follows: with I β a function of the masses of the particles running in the loop which contains three contributions corresponding to the three diagrams in Fig. 5. Thus, we express I β as follows: , (A.7) , (A.8) where we have passed to Euclidean momenta and integrated over the angular variables. Adding the three contributions we simplify the expression for I β as follows: . (A.10) We have checked that we obtain the same result by using the equivalence theorem where the external W -bosons are replaced by the corresponding Goldstone bosons.
In the limit m H = m A = m χ ++ = m χ + and m S m A we obtain I β ∼ 1/4 while if all masses are equal we get I β = 1/24. If we fix sin(α) m A can be obtained from m H and m S using Eq. (13b) while m χ + can be written in terms of m χ ++ and m A using Eq. (9). Thus, I β can be written as a function of sin(α), m χ ++ , m H and m S only. In Fig. 6 we present results for some representative values of the masses (we fix sin(α) = 0.08 and give I β as a function of m S for different values of m H = m χ ++ ).

Appendix B Details of the calculation of the neutrino masses
We define the Majorana mass matrix for the neutrinos as follows: Our parametrization for the elements of the neutrino mass matrix have been displayed in Eq. (24) which, in terms of the physical parameters, can be rewritten as In the unitary gauge there are four diagrams contributing to the neutrino masses as displayed in Fig. 7. As explained in Appendix A, each diagram will be finite when we add together the contributions from H, S and A. Note that the two diagrams in the last row of Fig. 7, after some relabeling of momenta, will give identical contributions. Taking this into account, we decompose I ν into three pieces as follows: κ ++ q P c 4M 4 W + M 2 W (q 2 1 + q 2 2 ) + (q 1 q 2 ) 2 (q 3 + q 1 + q 2 ) 2 + m 2 χ ++ , (B.4b) with, P c = 1 q 2 1 (q 2 1 + M 2 W )q 2 2 (q 2 2 + M 2 W ) (q 1 + q 2 ) 2 + m 2 κ ++ (q 2 3 + m 2 H )(q 2 3 + m 2 S )(q 2 3 + m 2 A ) , To evaluate the integrals in Eq. (B.4) we express the Euclidean four-momenta in the four dimensional spherical polar coordinates as follows: q i = q i (cos ψ i , sin ψ i cos θ i , sin ψ i sin θ i cos φ i , sin ψ i sin θ i sin φ i ) , (B.6) where, for brevity, we have used q i to denote both the four Euclidean vector and its modulus. With this, the differential under the integral can be expressed as: Without any loss of generality we can orient our 1-axis in the direction of q 3 and express the momenta as follows: q 3 = q 3 (1, 0, 0, 0) , q 2 = q 2 (cos ψ 2 , sin ψ 2 , 0, 0) , q 1 = q 1 (cos ψ 1 , sin ψ 1 cos θ 1 , sin ψ 1 sin θ 1 , 0) . (B.8) In this way, the integrands in Eq. (B.4) will not depend on the angles φ 1 , φ 2 , θ 2 , φ 3 , θ 3 , ψ 3 and they can be integrated out very easily. After this, the remaining six parameter integrals can be computed numerically (we have used Mathematica along with the Cuba package for this purpose). We have also checked numerically that, in the limit g → 0 and small mixing, our unitary gauge calculation agrees with the calculation discussed in Sec. 4, which includes only diagrams with scalar exchanges.
In Fig. 8 we give I ν as a function of m κ for different values of the other parameters. As in Sec. A we use Eq. (13b) and Eq. (9), fix sin(α) = 0.08 and take m H = m χ ++ ).

Figure 8:
The neutrino mass integral, I ν , as a function of m κ ++ for some representative values of the other parameters. We fix sin(α) = 0.08, use Eq. (13b) and Eq. (9) and take m H = m χ ++ .