Zooming in on eV-MeV Scale Sterile Neutrinos in light of Neutrinoless Double Beta Decay

The existence of light sterile neutrinos, as predicted in several models, can help to explain a number of observations starting from dark mater to recent anomalies in short baseline experiments. In this paper we consider two models - Left-Right Symmetric Zee model and Extended Seesaw model, that can naturally accommodate the presence of light sterile neutrinos in the eV to MeV mass scale. We perform a detailed study on the neutrinoless double beta decay process which receives major contributions from diagrams involving these light sterile neutrinos. Considering a number of theoretical and experimental constraints, including light neutrino masses and mixings, unitarity of the mixing matrix etc., we compare our predicted values of the half-life of neutrinoless double beta decay with the experimental limits. This can put significant constraints on the neutrino mass, active-sterile neutrino mixing and several other important parameters in these models.


Introduction
The Standard Model (SM) of particle physics, despite its major successes, is unable to explain the observed light neutrino mass splittings and their mixings, which provides a strong motivation to invoke beyond the Standard Model (BSM) physics.The two observed neutrino mass splittings are ∆m 2  12 ∼ 10 −5 eV 2 , |∆m 2 13 | ∼ 10 −3 eV 2 while the best-fit values of the neutrino mixing angles are θ 12 ∼ 34 • , θ 23 ∼ 48 • and θ 13 ∼ 8 • [1].Though neutrinos are massless in SM, a number of BSM theories have been proposed that successfully explain neutrino masses and mixings.One of the most appealing frameworks to generate Majorana masses of light neutrinos is via seesaw, where the dimension-5 lepton number violating operator generates the mass term after electroweak symmetry breaking [2][3][4][5][6][7][8].The type-I seesaw serves as the most economical framework, as the model in addition to the SM particles is minimally extended by gauge singlet right-handed neutrinos.Another popular class of mechanism is the radiative mass generation [9][10][11][12][13], where neutrino mass is generated via loop effect.In this work, we have considered a variation of the type-I seesaw model referred as Extended Seesaw model [14,15] and a left-right symmetric extension of radiative neutrino mass model [16][17][18][19].
The type-I seesaw model is the most economical, as the SM particle content is expanded with at least two heavy gauge singlet right-handed neutrinos which participate in light neutrino mass generation via seesaw mechanism.However, the drawback of this simplest model is that the mixing of these right-handed neutrinos with SM neutrinos are tightly constrained by eV light neutrino mass constraint, making the detection prospect of these right-handed neutrinos at experiments challenging.In Extended Seesaw, as the name suggests, more singlet neutrinos with large mixings are introduced with the possibility that some of them remain light and can be detected in experiments.The other popular mechanism for neutrino mass generation is through loop-induced processes.One of the simplest examples of this process is realized in the Zee model where the introduction of a doublet scalar and a charged singlet scalar can generate neutrino masses at the one-loop level.Though the simplest form of Zee model [20] cannot satisfy neutrino oscillation data [21][22][23], its left-right symmetric extension however is consistent with experimental observations [18,24].Here the Majorana mass of the left-handed and right-handed neutrinos are generated at the one-loop level while the Dirac mass term arises at the tree level from the Yukawa interactions.Finally the light neutrino masses are obtained by a type-1 seesaw mechanism but with the exception that the right-handed neutrinos can be light and hence offer better detection prospects.Both of the above mentioned models can accommodate light right-handed neutrinos with masses ranging in the eV to MeV scale.An eV scale sterile neutrino is well motivated, as this can explain LSND anomaly [25][26][27].Recently this anomaly has also been favoured by the MiniBooNE collaboration [28] but at the same time the data has been disfavoured by the KARMEN [29] and MINOS [30] observations.These issues may be finally tackled by the upcoming DUNE experiment [31,32].Further hints regarding the presence of an eV scale sterile neutrino comes from the reactor anti-neutrino anomaly (RAA) [33,34] and the Gallium anomaly [35,36].A keV scale sterile neutrino can be an excellent candidate for warm dark matter.Several disagreements between the cosmological observations and the N-body simulations of structure formations can be solved by introducing a keV scale warm dark matter candidate [37].The presence of an MeV sterile neutrino, on the other hand, can produce several observable astrophysical signals, such as, its effect on the Cosmic Microwave Background spectrum [38] and by producing X-ray photons which may be observable in satellite based X-ray experiments.
If the right-handed neutrinos are Majorana particles, they can give rise to additional contributions to the neutrinoless double beta decay (0νββ) process.The 0νββ process is the transition (A, Z) → (A, Z + 2) + 2e − with no neutrino being emitted [39][40][41].The process is lepton number violating [42,43].Depending on the mixing of the right-handed neutrinos with active neutrinos in type-I/Extended Seesaw, or the interaction of these right-handed neutrinos with the right-handed gauge boson in left-right symmetric extension, these right-handed neutrino states may give significant contributions in 0νββ process compared to the three SM neutrino contributions, and thus opening up the scope of detection of these Majorana neutrinos via 0νββ process.A number of experiments have searched for this process and the non-observation of signal has given bound on the half-life T 0ν 1/2 of 0νββ [44][45][46][47].The limit obtained on T 0ν 1/2 for 76 Ge is T 0ν 1/2 > 8.0 × 10 25 year from GERDA-II [44]; whereas, at 90% C.L. the KamLAND-Zen experiment has set a more stringent lower limit on the half-life of 136 Xe isotope as T 0ν 1/2 > 1.07 × 10 26 year [45].
In the present work, we consider Left-Right Symmetric Zee model (LRS Zee) and Extended Seesaw model that naturally accommodate light scale sterile neutrinos with masses ∼ eV to MeV.Our main goal is to study the 0νββ phenomenology for these two models.The left-right symmetric extension of the Zee model [18,24] presents a unique scenario where the model can be tested at the collider experiments as well as the neutrino experiments.It may have observable signals at the hadron [18,24] and lepton colliders [19] and most notably can be accessed at the very early stage run of the upcoming e + e − colliders.In addition to satisfying all the neutrino mass and mixing constraints, this model can also give rise to several new 0νββ process which can significantly enhance the decay rate.This results in a marked decrease in the half-life of the 0νββ decay process in this model.We study the variation of T 0ν 1/2 with respect to different model parameters and identify three of them which are most significant.These three parameters are the lightest neutrino mass, the Dirac CP phase and the mixing angle between the left and right gauge bosons.By varying these parameters we identify the regions which can be ruled out from the experimental limits on the half-life of 0νββ process.In the case of extended seesaw mechanism, we first give an approximate analysis considering only one generation right-handed neutrino and one generation active neutrino.Subsequently, we present a realistic analysis of the half-life with three generation active neutrinos and six generation right-handed neutrinos.We have considered all the constraints arising both from theory and experiments.For the active neutrinos, we have considered bounds on the mass square differences, three mixing angles in agreement with neutrino oscillation data [1] and the limit on the sum of the masses of active neutrinos which comes from Planck satelite experiment [48].We have ensured a mass hierarchy among these active and right-handed neutrinos to validate the seesaw approximation for this model, as well as have considered constraints from non-unitarity [49].We have calculated the 0νββ decay contribution considering all the required model parameters which pass all the aforementioned constraints and have checked if the predicted contribution satisfies the corresponding experimental limits [45].
The paper is organised in the following way.In Section 2 we present a detailed study of the 0νββ process in the LRS Zee model.This is followed by Section 3, where we give a detailed description of extended seesaw scenario and the analysis of the model with respect to many theoretical and experimental aspects, e.g., neutrinuo oscillation data, 0νββ decay, unitarity and others.Finally, we present our conclusions in Section 4.

Left-Right Symmetric Zee Model and Analysis
The Zee model [20] is one of the simplest extensions of the Standard Model (SM) which can explain the origin of neutrino mass.By extending the SM framework with an extra doublet and a charged singlet scalar, neutrino masses can be successfully generated at the one-loop level.This simplest form of the Zee model, though, is found to be incompatible with the neutrino experimental data [21][22][23]50] and one needs to extend it further in order to get a viable scenario to explain all the neutrino oscillation constraints.The left-right symmetric extension of the Zee model [18,24] provides an alternate model framework which can easily explain the neutrino oscillation data as well as provide interesting flavor violating signals and unique collider signatures.
The pair production and decay of the singly charged Higgs boson can produce final states with two charged leptons (with either same or different flavors) and missing transverse energy.This process has been studied in the context of the Large Hadron Collider (LHC) in Ref. [18,24] and for the International Linear Collider (ILC) and the Compact Linear Collider (CLIC) in Ref. [19].The charged singlet scalar pair-production cross-section at hadron colliders is quite small and is dominated by the photon mediated process, while it can become significantly larger in the electron-positron collider due to the right-handed neutrino mediated t-channel diagram.Thus the ILC and CLIC experiments may be able to observe such a particle with very low integrated luminosity of only 1-3 fb −1 , see [19] for details.Below we present a brief discussion on the model.

Model
Other than the three SM neutrinos, the neutrino sector in this model also contains light right-handed neutrinos with masses ranging from an eV to around 100 eV.Their masses depend on the coupling of righthanded neutrinos with the charged singlet scalar.The Yukawa Lagrangian in this model is given as [18,19,24] where Φ and δ are the bidoublet and charged singlet scalars with Y and λ being their respective Yukawa couplings.The Majorana masses of all the neutrinos are generated at the one-loop level and as a result, they remain quite light.In order to understand the neutrino sector, we also need to understand the scalar sector of the model.

Charged scalar spectrum:
The minimal Higgs sector in this model consists of a bidoublet, two doublets and a charged singlet field given as The SU (2) R × U (1) B−L symmetry is broken down to U (1) Y as the neutral component of the right-handed doublet H R acquires a non-zero vacuum expectation value (VEV).The bidoublet Φ is required to generate the quark and charged lepton masses and Cabibbo-Kobayashi-Maskawa (CKM) mixing angles.The neutrino masses, similar to the Zee model, are generated at the one-loop level due to the presence of the charged singlet scalar field.The charged Higgs bosons play an important role in the generation of the one-loop neutrino masses as their masses and mixings become important parameters in the expression for the induced neutrino Majorana masses.So before we study the neutrino mass generation mechanism it is important to define the mass basis for the charged Higgs bosons.There are in total five charged Higgs states which mix to give five mass eigenstates through the rotation where V is the 5×5 charged Higgs mixing matrix.There are three physical Higgs bosons H + 1 , H + 2 , H + 3 which will contribute to the neutrino masses and two Goldstone states G + 1 , G + 2 which are eaten up by the W R and W boson as their longitudinal degrees of freedom.In this model, the φ − 1 * and the H + R become the Goldstone bosons while the other three charged states mix to form the physical Higgs bosons.Constraints from flavor violating process further require the mass of the second bidoublet scalar φ + 2 to be heavier than around 15 TeV.Thus the largest off-diagonal contributions to V comes from the mixing of the charged singlet δ + with the left-handed charged Higgs boson H + L .These elements will play an important role in the neutrino mass generation as well.Neutrino mass and mixings: The neutrino sector consists of three left-handed and three right-handed neutrinos.The absence of triplet scalars in the model prevents us from writing a Majorana mass term for the neutrinos.All the neutrino Majorana masses here are generated at the one-loop level and hence are quite small.The lightest right-handed neutrino mass ranges from a few eV to a few hundred eV and the other right-handed neutrinos also remain lighter than an MeV.The Dirac masses are thus required to be quite small as well, so as to satisfy the experimentally observed neutrino masses and mixings1 .The one-loop Feynman diagram for the generation of neutrino Majorana masses is given in Fig. 1.The corresponding expressions for the neutrino Majorana masses in this case are given as [24]: In the above, λ L/R αβ = λ αβ L/R − λ βα L/R , m e β and M hi are the charged lepton and Higgs boson masses and V ij are the charged Higgs boson mixings given in Eq. 3.For our calculations we will consider the case with λ L = 0 as was discussed in [19].This will thus give us a 6 × 6 neutrino mass matrix given as where We thus have a scenario that is very similar to the type-I seesaw mechanism, i.e., the light and heavy neutrino mass matrix after persuing a block-diagonalization becomes, The neutrino rotation matrix, taking it from flavor to mass eigenstates, will be a 6 × 6 matrix which we can write as such that Here M ν = diag(m 1 , m 2 , m 3 ) and M N = diag(M 1 , M 2 , M 3 ) are diagonal matrices consisting of the light and heavy neutrino masses, respectively.Fig. 2 shows a plot of the Majorana masses of the right-handed neutrinos as generated at the one-loop level in this model.The plot shows the variation of the eigenvalues of the right-handed neutrino Majorana mass matrix (elements of the matrix given in Eq. 4) as a function of its coupling with the charged singlet scalar.As λ R is varied from 0.1 to 3, the lightest right handed neutrino mass M N1 varies from 3 eV to 80 eV, while M N2,3 remain in the sub-MeV scale.In generating the plot, we considered the lightest charged Higgs boson has a mass 473 GeV, which primarily consists of charged singlet δ.Below we discuss the contribution of the right-handed neutrinos in 0νββ decay.

Diagrams and amplitudes of 0νββ transition
Contrary to most seesaw models which contain the right-handed neutrinos of TeV scale mass, the LRS Zee model naturally accommodates eV-MeV scale right-handed neutrinos, as has already been discussed in the previous section.Diagrams involving the right-handed neutrinos can thus significantly contribute to the 0νββ processes and this model gives us an excellent framework to study these effects.
The Feynman diagrams of all the possible contributions are presented in Fig 3 .For each diagram we write its amplitude and identify the dimensionless parameter η i that will be used in the computation of the half life (T 0ν 1/2 ) of the 0νββ process.In the subsequent discussion, we refer the mass eigenstates of SM neutrinos as 'light', and the right-handed neutrinos as 'heavy', as the right-handed neutrino states are heavier than the SM neutrinos.
• Light neutrino diagram: The diagram (a) corresponds to the light neutrino contribution.Its amplitude is given as where G F is the Fermi constant, p is the momentum transfer at the leptonic vertex and i = 1, 2, 3 corresponds to the light neutrino mass eigenstates.The corresponding η obtained in this case is given as • Heavy neutrino diagrams: Diagram (b) corresponds to the heavy neutrino contribution.The heavy neutrinos in this model are composed of the right-handed neutrinos but unlike other Left-Right models, they are quite light in this case with masses in the eV to MeV range.The Feynman amplitude and the corresponding η from this diagram is given as • λ diagrams: Diagrams (e) & (f) represent the processes mediated by the W L − W R exchange.The Feynman amplitudes from each diagram can be easily combined to give us a final expression which is and the expression for the η parameter is • η diagrams: Diagrams (g) & (h) are due to the W L − W R mixing in this model and depend on the The Feynman amplitude combining the two diagrams can be written as and the corresponding η parameter is The half-life for the 0νββ process after combining the contributions from all these diagrams is then given as [51-54] where G 0ν is the phase space factor; M 0ν ν , M 0ν λ and M 0ν η are the nuclear matrix elements.Now that we have the expression for the half-life for the 0νββ processes, let us discuss some of the features of this framework which will help us understand the relative contribution arising from each of these diagrams.The right-handed neutrino masses being at the eV to MeV scale contribute significantly to these processes here and hence the diagrams involving N R become quite important.The relative contributions from the diagrams are also highly dependent on the light-heavy neutrino mixings (S, T) with the λ and η diagrams becoming significant as this mixing increases.The gauge boson (W L − W R ) mixing is another important factor in these diagrams and its value can determine which diagram gives significant contribution to the 0νββ decay process.Finally since the W R boson mass is required to be quite large from experimental constraints [55][56][57][58][59][60][61][62][63][64][65][66][67], we have chosen it to be 5.5 TeV here and this results in a large suppression for all the diagrams with amplitudes involving (M W L /M W R ) 4 term.
The neutrino parameters in this model depend significantly on the masses and mixings of the charged scalars as can be seen quite clearly from Eq. 4. A close inspection of this equation also shows that the factor V 5i which is the mixing between the charged singlet and other charged Higgs states is quite important for the neutrino masses.As discussed earlier, the charged singlet Higgs δ + can only have significant mixing with the left-handed charged scalar H + L .We can thus approximately write where H + 1 and H + 2 are the lightest and next-to-lightest charged Higgs bosons respectively with θ being the mixing angle.Clearly two extreme cases appear here 1. maximal mixing with θ = π/4, denoted as H max , 2. minimal mixing with θ = 0, denoted as H min .

Results
To analyse the half-life of 0νββ process for Germanium ( 76 Ge) and Xenon ( 136 Xe) nucleus, we consider two cases of maximal and minimal mixing as described in the previous section.We consider a normal mass ordering among the SM neutrinos, and use the Casas-Ibarra [68] parametrization to fit the latest neutrino oscillation data [1].As the right-handed neutrino masses are quite small here, the mixing between them and the left-handed neutrinos, represented by the S and T matrices in Eq. 7, can become quite significant 2 .For a fixed choice of the right-handed neutrino masses, the light-heavy neutrino (ν L − N R ) mixing depends largely on the light neutrino masses.As the SM neutrino masses increase, the mass difference between the light and heavy neutrino states become smaller resulting in a larger mixing angle.Thus the lightest neutrino mass m ν1 is an important parameter for our analysis.
The other parameters which play a significant role in determining the value of T 0ν 1/2 are the W L −W R mixing angle (θ LR ) and the Dirac CP phase (δ CP ) of the neutrino Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix.The contribution from the η diagram, directly proportional to tan θ LR , can become substantial depending on this mixing angle.The value of δ CP , though directly does not appear in any of the expressions, 2 These terms still remain orders of magnitude smaller than the U and V matrices.
determines the neutrino parameters obtained from the Cassas-Ibarra parametrization.This has a significant consequence on the calculated value of T 0ν 1/2 .The nuclear matrix elements (NMEs) for 76 Ge and 136 Xe, which we adopt from [44,45] are equally important for the evaluation of T 0ν 1/2 .We consider two cases, one with the maximum and another with the minimum values of the NMEs, and evaluate the half-life.For  For each of these cases we vary m ν1 , θ LR and δ CP to obtain the predicted value of the half-life of 0νββ decay process.Fig. 4 shows the variation of T 0ν 1/2 for 76 Ge nucleus as a function of the lightest neutrino mass  The shaded region correspond to T 0ν 1/2 < 8.0 × 10 25 years for left panel, and 1.07 × 10 26 years for right panel and disallowed by GERDA [44], and KamLAND-Zen [45], respectively.m ν1 for a fixed value of θ LR and δ CP .The values of all other PMNS matrix elements were fixed to their central values and the λ R matrix was chosen such that the right-handed neutrino masses were 7.92 eV, 3.54 keV and 3.55 keV respectively.As can be seen here, the half-life decreases quite drastically as the lightest neutrino mass increases.This is because the light-heavy neutrino mixing increases as discussed earlier and as a result the η 3 , η λ and η η contributions become dominant.As for this figure, we consider a large value of θ LR , therefore η η always dominate.We find that the canonical light neutrino contribution η 1 is rather subdominant in this figure.Fig. 5 shows the variation of T 0ν 1/2 with the left-right charged gauge boson mixing θ LR .As θ LR increases the decay half-life falls drastically for a value θ LR 10 −6 .This is the point at which the η η term, which is proportional to tan θ LR starts dominating over the other terms resulting in a steep decrease of the half-life as expected.For smaller value of θ LR , the dominant contribution arises from η 1 , η 3 , η λ , which are independent of the left-right mixing.Since the CP violating phase δ CP is another crucial parameter in our analysis, we show the variation of half-life with respect to δ CP .Fig. 6 gives the change of T 0ν 1/2 for 76 Ge nucleus as a function of δ CP .A close inspection of the numbers we obtained shows that the variation of half-life mirrors the variation in sin δ CP which directly determines the values of the neutrino parameters obtained in our calculations.This is to note that, in all these figures, the scenario H min − Ge max gives the strongest constraint.This can be understood from the expression of T 0ν 1/2 as given in Eq. 18.As can be seen, the calculated half-life is smaller (leading to a more constrained scenario) for larger values of the amplitudes and NMEs.Therefore, naturally the maximum values of 76 Ge NMEs leads to a more constrained scenario.
The charged Higgs mixing on the other hand plays an indirect but significant role in determining the values obtained for the Feynman diagram amplitudes.As was discussed earlier, the Feynman amplitudes corresponding to η 2 , and η 4 are negligible due to the (M W L /M W R ) 4 suppression.So the dominant contribution always arises from any one of the η 1 , η λ or η η terms, while we find η 3 contribution is slightly smaller than these above mentioned contributions.A smaller charged Higgs mixing will invariably lead to a lighter Majorana mass for the RH neutrinos, which has a two-fold effect on the neutrino sector.Firstly, lighter RH neutrinos will results in relatively heavier active neutrinos since the active neutrino mass is obtained by seesaw mechanism in our case.This will boost the η 1 amplitude resulting in a smaller value of T 1/2 .Secondly, a heavier active neutrino will result in a larger light-heavy mixing as discussed earlier.This again helps boost the η λ and η η amplitudes further lowering the calculated half-life of 0νββ process, leading to a tight constraint on the parameter.The plots obtained for the 136 Xe nucleus are very similar in nature to the ones for 76 Ge and the most constrained scenario is again the H min − Xe max case.This warrants for a more detailed study of this particular case for a better understanding, which we discuss below.We present the results for both 76 Ge and 136 Xe nucleus in the ensuing discussion of the most constrained scenario for each, i.e., largest values for the NMEs and minimal mixing of the Higgs sector H min − Ge max and H min − Xe max .Fig. 7 shows the variation of T 0ν 1/2 for 76 Ge as a function of the lightest neutrino mass for several fixed values of θ LR and δ CP .It is quite clear that the value of T 0ν 1/2 decreases as the lightest neutrino mass and/or value of θ LR increases.The variation with m ν1 is quite evident from our earlier discussion.For the set of plots with θ LR = 10 −8 , the initial variation with m ν1 is quite modest till m ν1 10 −7 eV.The major contribution here comes from the η 1 term with the half-life slowly decreasing with an increase in the lightest neutrino mass.At larger values of m ν1 , the light-heavy neutrino mixing increases significantly and the dominant contribution comes from the η λ term.The effect of increasing θ LR can be understood as an artefact of an increase in the η η term which starts contributing quite significantly at θ LR 10 −6 .Similar characteristics and dependence can be observed in Fig. 8 where we have plotted the variation of T 0ν 1/2 as a function of θ LR for fixed values of m ν1 and δ CP .As can be seen here, the value of T 0ν 1/2 remain almost constant in the initial region of θ LR 10 −6 .In this region, the dominant contribution comes from η 1 and η λ and since the lightest neutrino mass is constant for each line, there is no variation in their value.In the region θ LR 10 −6 , the η η term starts dominating and we see a sharp decrease in the value of T 0ν 1/2 in this region.
Finally, the variation of T 0ν 1/2 with δ CP is given in Fig. 9.This plot is quite interesting as it clearly shows the contribution of different η terms in different regions of parameter space.The line corresponding to θ LR = 10 −8 and m ν1 = 10 −9 eV has a dominant contribution from η 1 term.As η 1 = 1 me i U 2 ei m i , its contribution in this case only depends on the matrix elements U ei since m e and m i remains constant.The U 13   element 3 is directly proportional to e −iδ CP and as a result one can approximately write As δ CP approaches 180°or 360°, the value of η 1 goes through a maxima while it becomes a minima at δ CP = 270°.The inverse of this behaviour is reflected in the T 0ν 1/2 plot.For the other lines in this figure, they receive dominant contributions from either η λ or η η .Here as U 13 decreases, the elements of S and T mixing matrices increase and hence their nature is opposite to the previous plot.The nature of the plots for 136 Xe nucleus is the same as the 76 Ge nucleus except the fact that the parameters are more tightly constrained owing to the more stringent experimental limit of the 0νββ half-life for the 136 Xe nucleus. 3The U 13 element will be the same as the (1,3) element of the PMNS matrix.Fig. 10 shows the upper limits on the mass of the lightest neutrino in the this model as a function of θ LR for fixed values of δ CP for both 76 Ge and 136 Xe nucleus.As we have already discussed that the most stringent bound on the parameter space is obtained for a δ CP of around 264°, this fact is also reflected from this figure.As expected from the previous discussion, the upper limit on the lightest neutrino mass becomes much stronger for larger values of θ LR and vice-versa.This occurs due to dominant η η contribution for a large θ LR , leading to a tighter constraint on the lightest neutrino mass m ν1 .Another observation is that the limits obtained for the 136 Xe nucleus is much stronger than the 76 Ge nucleus.

Extended Seesaw Model and Analysis
This is another extension of SM, where the model contains light and heavy sterile neutrinos 4 , which can give large contribution in 0νββ process.Several studies [42,43,51,53,[69][70][71][72][73][74][75][76][77][78][79][80][81][82][83] have focused on sterile neutrinos with mass > 100 MeV, and large contribution in the LNV process.Here, instead we consider some of the sterile neutrino states in the < 100 MeV mass range.We investigate the allowed model parameters, which satisfy a number of theoretical and experimental constraints.In doing so, we first consider a simplistic one generation scenario with one active neutrino, one light sterile neutrino S L , and another heavy sterile neutrino N L .Subsequently we extend our analysis with realistic three-generation case where the neutrino sector comprises of three active neutrinos along with the additional six sterile (three S L and three N L ) which are neutral under the SM gauge group.Below, we first review the model, and then discuss the contribution in 0νββ process.

Model
The neutral lepton sector of the model contains three generation of SM neutrinos ν L along with additional sterile neutrino states S L and N L .The mass term of the neutrinos have the following form, We denote the neutral lepton mass matrix as M n , and hence We choose to work in a basis where the Majorana mass matrix M R of N L sterile is real.The term containing µ denotes the Majorana mass of the heavy neutrino state S L with µ being a complex symmetric matrix.The Dirac mass matrix M D represents the mixing between the SM neutrino states ν L and the heavy sterile neutrino states N L ; whereas, M S represents the mixing between the two sterile states S L and N L .Throughout our analysis we consider the matrices M R and M S are invertible.We also assume, that, the different sub-matrices follow the hierarchy, M R > M S > M D µ.For seesaw approximation to be valid, the mixing matrices should satisfy [15] for details.
Contrary to the inverse seesaw [84][85][86][87][88][89][90][91][92][93][94][95], Extended Seesaw model contains both the heavy and small lepton number violation scales M R and µ respectively.The SM neutrino masses strongly depend on the small lepton number violating scale µ and hence in the µ → 0 limit, the ν L states become massless.The heavy Majorana neutrino contribution in 0νββ decay can be sizeable, even in the µ → 0 limit.Hence, the contributions of the SM neutrinos and the heavy Majorana neutrinos in 0νββ process are completely decoupled from each other.Contribution from heavy sterile neutrinos for this model has been discussed in [15].
The neutrino mass matrix M n can be diagonalized by a unitary transformation, where U as an expansion with order parameter M D /M S has the following form [15]: In the above, W µ , W S and W N are the three unitary matrices that diagonalize the block diagonal matrices The matrix m ν represents the light neutrino mass matrix, and m s and m n represent the heavy neutrino mass matrices.The hierarchy among the sub matrices ensures that m n and its eigenvalues give the heaviest sterile neutrinos in this model.The other sterile neutrinos that originates from the diagonalization of m s can be relatively lighter, but they certainly should be heavier than the three active neutrinos m ν < m s < m n .In the subsequent sections, we explore the possibility of that the sterile states from m s are in the eV to MeV range, while the remaining sterile neutrino states m n are more than GeV.Before presenting a detailed analysis on 0νββ, we first consider additional constraints coming from light neutrino mass measurement, non-unitarity and others.

Constraints
Before delving into the analysis, we present a short descriptions of all constariants that has been applied in this model.

(a) Theoretical Constraints:
Hierarchy: The different sub-matrices of Eq. 20 should satisfy the hierarchy , m s > µ (from Eq. 24), for one generation this will be M 2 S /M R > µ.This limit also defines the region where seesaw approximation is valid [15].

Unitarity:
The mass matrix being symmetric, the diagonalization matrix given in Eq. 22 should be orthogonal or unitary, i.e., we should have the relation U † U = UU † = I; but working with the seesaw approximation up to a certain order and also having low scale sterile, the elements of U = UU † will not be identity matrix rather those elements will be I ± δ, where δ is the tolerance of every single elements of U to get a viable parameter space for lightest sterile 5 in this model.So, to zero in on the allowed parameter space of eV to MeV sterile in this model, we have to constrain the parameter space, setting some cut-off values on the both diagonal and non-diagonal elements of U .In short, we allow some error bar on the diagonal elements on U around unity and for non-diagonal elements the required error bar will be around zero.Depending on the choice of parameter space, the error bar may differ for diagonal and non-diagonal elements.We have generally taken the the maximum constraints on the deviation which provides us the desired allowed parameter space.
(b) Experimental Constraint: Mass of active neutrino: We consider the constraint on the sum of active neutrino masses from Planck cosmological data [48], i.e., at 95% C.L. the sum on the masses of active neutrinos will be less than m ν < 0.194 eV.In the analysis of one-generation, this bound simply manifests as the upper bound on the mass of single active neutrino.We implement the constraints on mixing angles and on the mass square differences among three active neutrinos from neutrino oscillation data in the three-generation case [1,96] in case of Normal hierarchy.
Constraints from 0νββ limit: The limit on the T 0ν 1/2 from the KamLAND-Zen [45] severely constrains the parameter space for eV/MeV sterile of this model, see Sec. 3.3.

Daya Bay experiment:
The Daya Bay reactor anti neutrino experiment shows a large exclusion region between 2 × 10 −4 < ∆m 2 s1 < 0.3 eV 2 as function of sin 2 2θ 1s [97] at 95% CL, where ∆m 2 s1 is the mass-square difference between extra sterile and electron-neutrino (ν e ) and θ 1s is the angle of active-sterile mixing.This result will further constrain the allowed parameter space for eV sterile.

0νββ decay: sterile neutrino contributions
In this section, we outline the contributions of sterile neutrinos having Majorana masses in 0νββ decay.
The half-life of 0νββ is written as [15,98] 1 where j represents the number of light neutrino states and the additional heavy neutrino states.The parameters µ j and Θ ej represent the masses of the neutrino states and the mixing with SM neutrinos respectively.In the above, Mν .The reference mass scales are considered as electron (m e ) and proton (m p ) masses, M ν and M N are the NMEs for exchange of light and heavy neutrinos respectively.The values of NME and phase space factor G 0ν have been taken from Ref. [99].Below, we classify the sterile neutrino contributions according to its mass scale.
Other than the contributions of the SM neutrinos, the sterile neutrino states S k and N k (k = 1, 2, 3; in our case) can also contribute in the 0νββ process.Evidently, we have two extra contributions apart from the SM one.
(a) The heaviest states N k have a mass range m n k > 100-200 MeV and give a contribution in 0νββ as where A N represents the amplitude of this process, and V eN k is the mixing of the N k states with the active neutrinos.Using V eN = M † D M −1 R W N , this can be simplified as (b) The other sterile neutrino states S k give contributions proportional to for the mass range m s k > 100-200 MeV, whereas when sterile mass is light.We use the compact expression for the amplitude, that also take into account The value of m n k in our analysis is 10 5 GeV for one-generation case giving rise to active-sterile mixing as M D /M R ∼ 10 −10 − 10 −7 , whereas m n k s are to be of the order of 50 to 500 GeV in the three-generation one having V 2 eN k as ∼ 10 −8 to 10 −7 .So, the sterile neutrinos N k being heavy contributes negligibly and hence we do not consider its contribution.The half-life of 0νββ is thus, where j represents the number index of three light neutrino states whereas k denotes the number index of comparatively light sterile states S L .This is to note that, the lighter S L states (eV-MeV) can have substantial active-sterile mixing To get the essence of all theoretical and experimental constraints properly, first we consider only one generation for all states, so the mixing matrix given in Eq. 23 will be of the order of 3 × 3 having M D , M R , M S , µ as simply numbers.

One Generation
In this section, we provide a detailed analysis and results of the allowed range of MeV/eV sterile neutrino.In Fig. 11, the plot in the left panel (Fig. 11a) shows the allowed region for MeV sterile and its contribution to 0νββ in M D −M S plane for one-generation scenario.The matrix M n in Eq.21, is 3×3 instead of being 9×9 in this case.The square boxes in the index box of this figure (Fig. 11a) shows the color of allowed regions in agreement with different constraints and the respective texts are written in the same color of the border of that region.The cyan colored region enclosed by the red-dashed curve in Fig. 11a shows the region allowed by GeV.The orange shaded region of (b) will be constrained further from reactor anti-neutrino Daya Bay experiment [97] that will be elaborated in Fig. 12.
the off-diagonal element U 13 , where U 13 is the (1,3) element of U † U with U being the diagonalization matrix (Eq.22).Here, we consider the constraint as |U 13 | < 10 −8 , i.e., U 13 is almost vanishing.The lower green region enclosed by blue solid line presents the region allowed by the diagonal element U 33 .The constraints due to other matrix elements of U coming from the condition of diagonalization matrix U being unitary are less stringent and are simply allowed by the final overlapped zone.In that final allowed region, the numerical value of off-diagonal element U 13 is ∼ 10 −9 and that of diagonal element U 33 is ∼ (1 + 10 −16 ).The region covered by pink-colored straight lines shows the mass of light neutrino range 0 < m ν < 0.194 eV.The brown rectangle region enclosed from below by the black dashed line (near M S = 32 GeV) represents the constraint m s = M 2 S /M R > µ marking the area where the seesaw approximation is valid.The extreme left almost vertical orange-colored region enclosed by the solid orange bow-type curve shows the region in agreement with the contribution of 0νββ, where T 0ν 1/2 > 1.07 × 10 26 yr [45].In obtaining this allowed parameter space, we have considered both the light neutrino and sterile neutrino contribution, see Eq. 31.The values of the NMEs that we have considered in this analysis are M ν = 2.29 and M N = 163.5 [99].
The overlapped region in Fig. 11a enclosed by dashed black straight line from below and solid blue line from above with 32 GeV < M S < 49 GeV and red-dashed line from the right is the final allowed range for MeV sterile in Extended Seesaw model, with the value of M D ≤ 0.00011 GeV for M S ∼ 32 GeV and with M D value up to ∼ 0.002356 GeV for M S = 49 GeV.The allowed mass range of m s is 10 MeV < m s < 24 MeV.The mass of the active neutrino in that region is m ν ≤ 10 −2 eV.
The plot in the right panel (Fig. 11b) of Fig. 11 shows the allowed region for eV sterile in the Extended Seesaw model and its contribution to 0νββ decay in M D − M S plane.The inclusion of sterile neutrinos whether being heavy or light has its effect on the unitarity of PMNS matrix [49].The PMNS matrix Since, in Extended Seesaw, we have eV sterile starting from ∼ 0.1 eV, the region can further be constrained from reactor anti-neutrino experiments, such as, Daya Bay.Evidently, the orange shaded region of Fig. 11b or a part of it, where we have sterile ∼ O(eV) can be probed again from such results.We represent the mentioned region of Fig. 11b  The region in Fig. 11b, below the green line and right side of the red-dashed line being completely ruled out from our model parameters (Fig. 11b) represents the white regions of Fig. 12.The dotted red and black solid line represent the Daya Bay experimental constraints on low-scale sterile in ∆m 2 s1 (eV 2 ) − θ 2 1s plane [97].The gray-colored region shows the zone that are not allowed by this experimental data.The overlapped region of this adjacent figure is ruled out from the experimental result.We can see that Daya Bay results exclude some region from the left side (giving constraint on M D ) but still allow all the values of sterile neutrino from 0.1 eV to 0.4 eV.The mass of active neutrino in the remaining allowed zone is m ν ≤ 10 −4 eV.
In passing, we would like to comment on the active-sterile mixing value which is getting constrained from Daya Bay data.The θ 2 1s is actually (M D /M S ) 2 .Also, from unitarity [49], we have Fig. 12 shows for some M S , M D values θ 2 1s < 10 −3 giving slightly more stringent bound on the mixing compared to that of coming from unitarity in our model set-up.

Three Generation
In Section 3.4.1,we have discussed different constraints from neutrino mass, half-life of 0νββ decay, unitarity and from the validation of seesaw approximation for one generation realisation of Extended Seesaw model.In this section, we are extending the analysis for three-generation case which is more realistic than the previous scenario.In addition to the bounds from 0νββ, other experimental and theoretical constraints, we also satisfy neutrino oscillation data.In particular, we consider that, • upper bound on the sum of all three active neutrinos is constrained from cosmology, i m νi < 0.194 eV at 2σ C.L. [48], • two mass squared differences 6.93 < ∆m 2 21 10 −5 eV 2 < 7.97 and 2.37 < ∆m 2 31 10 −3 eV 2 < 2.63 vary in the 3σ range [1], • 3σ range [1] of the three mixing angles 30 In the present set-up, contribution to the 0νββ can come from light active neutrinos (m νi , i = 1, 2, 3), additional eV to MeV scale sterile neutrinos (m si , i = 1, 2, 3) and heavy GeV scale neutrinos (m ni , i = 1, 2, 3).As shown in Sec.3.2, the contribution of heaviest sterile neutrinos to 0νββ is evidently suppressed.Here, the neutrino mass matrix is 9 × 9, but we are working in the seesaw approximation regime which gives three sets of 3 × 3 matrices namely m ν , m s and m n .After diagonalization of each 3 block of Eq. 24 individually, we check the unitarity constraints as described in Sec.3.2 and in Sec.3.4.1.The matrices which diagonalize each blocked-matrices combine to form a 9 × 9 matrix (Eq.23) and we impose constraints on the unit matrix (U , see Sec. 3.4.1)with the absolute variation of each elements by ±10 −2 .This ±10 −2 variation manifestly impose constrains on the ratio M D /M S .We have checked that the error bar is consistent with the experimental bound [49] that arises due to the non-unitarity effect for eV-keV scale sterile neutrino.A detailed description of different conditions provided in Sec.3.2 to constraint the parameter space have been thoroughly followed in the present analysis.
As discussed earlier, after using the seesaw approximation, we get three different 3 × 3 matrices, which are m n , m s , m ν .Among them one corresponds to the mass matrix for the three active neutrinos (denoted by m ν ), other two correspond to the mass matrices for the three relatively light sterile neutrinos (denoted by m s ) and the three heavy neutrinos which are in GeV scale (denoted by m n ).From Eq. 24, we can see the expressions of m ν , m s and m n depend on the matrices M D , M R , M S and µ whose elements are the free input parameters in the extended seesaw scenario.We choose the model parameters in our framework in a way so that we can accommodate eV to MeV scale sterile neutrino.Before proceeding, we consider few assumptions which include M D , M S , M R as the real diagonal matrices and µ as the complex symmetric matrix i.e. µ T = µ and µ * = µ.In this work, we have focused on the normal hierarchy of the neutrino masses as illustrative example.In order to satisfy neutrino oscillation constraints and to obtain sterile neutrinos in eV to MeV scale, we have varied the model parameters in the following range (in GeV), where i, j can vary from 1 to 3. Below we show the allowed model parameters as well as correlation between different observables for this model as scatter plots.
In the left panel of Fig. 13, we have shown the allowed region in m νi (eV) − m s1 (GeV) plane after satisfying all data 6 , where m s1 is the physical mass of the lightest sterile neutrino state.In the figure, green dots show the range allowed by neutrino oscillation data (N.O.D), blue rhombus points represent the range allowed by 0νββ and unitarity along with N.O.D. Finally the red points exhibit the range that is being further constrained by m si > µ d i .Here, m si are the physical masses of the sterile state S L and µ d i are the eigenvalues of µ and i =1,2,3.In the present work, the model parameters are less constrained from the 0νββ decay bound than the unitarity and m si > µ d i bounds.The model parameters range considered in this work give us eV to MeV scale sterile neutrino as seen by the range of m s1 -axis.One interesting thing to note here is that m s1 ≥ 10 −6 GeV is disallowed when we consider both the constraints, unitarity and 0νββ.This is mainly due to unitarity bound since this bound mostly depends on the ratio of M D /M S .Therefore, when (M D /M S ) 2 < 10 −2 , those points satisfy the unitarity constraints which is ±10 −2 variation around the unit matrix (see Sec. 3.2).The disallowed points correspond to higher ratio, i.e., (M D /M S ) 2 > 10 −2 and those points represent lower values of the elements of µ matrix in order to satisfy the N.O.D which is not covered in Eq. 33.This also implies that the elements of M D and M S are of the same order for the disallowed points and more likely to have higher M S values.Finally the red points are obtained when we impose the In the left and right panels of Fig. 16 M S11 plane respectively.Here, the superscript I denotes the imaginary part.In the left panel we can see that most of the points which satisfy oscillation data are below the yellow line which corresponds to As we have discussed, the unitarity bound (variation of ±10 −2 around unit matrix) mostly depends on the ratio M Dii M Sii (i = 1, 2, 3) and we can roughly say if M Dii M Sii < 10 −1 , then it can pass the unitarity bounds.One interesting thing to note here is that there exist a sharp correlation among the M D22 and M S22 parameters (which is valid for other elements of M D and M S matrices also).This is because, approximately (M D /M S ) 2 µ ∼ 10 −11 GeV and we have taken µ < 10 −8 GeV, so M D22 and M S22 can not take arbitrary values which corresponds to significant difference in their magnitudes, otherwise neutrino mass data will not be satisfied.On the other hand in the right panel of Fig. 16, we can see that after imposing all the constraints we get the points which are more prone to have higher values of µ I 23 .This is because when we impose the unitarity constraint which corresponds to M D ii M S ii < 10 −1 (i = 1, 2, 3) and from the order of magnitude estimation of neutrino mass we obtain µ 10 −9 .This is clearly reflected in the right panel of the figure because the points are more dense in that region where µ I

Conclusion
In this work, we consider two theory frameworks with sterile neutrinos -a) Left-Right Symmetric Zee, and b) Extended Seesaw model, which successfully explain the light neutrino masses and their mixings.Both of the models can accommodate sterile neutrinos with their masses being free parameters varying over a wide range.We particularly focus on relatively lighter mass range ∼ eV to MeV, and explore the contribution to 0νββ process.The Left-Right Symmetric Zee model represents a scenario where the masses of the sterile neutrinos are generated at the one loop level.They are directly dependent on the right-handed Yukawa coupling λ R and the masses and mixings of the charged Higgs bosons.For large values of λ R (close to the perturbative limit) the sterile neutrino masses always remain well below the MeV scale.This presents a unique scenario where the three light sterile neutrinos can have significant contribution to 0νββ process.We find that the half-life of this process crucially depends on three parameters -lightest neutrino mass m ν1 , Dirac CP phase δ CP , and the W L − W R mixing angle θ LR .In our analysis, we consider the cases with maximal and minimal mixings among the charged Higgs bosons and also consider both the upper and lower values of the NMEs for 76 Ge and 136 Xe nuclei.The scenario with minimal mixing of the Higgs bosons and maximum values of the NMEs produces the most stringent bound on the model.The calculated half-life for both 76 Ge and 136 Xe nuclei decreases drastically with an increase in m ν1 or θ LR .This is due to the dominant contributions coming from the λ and η diagrams as m ν1 and/or θ LR are increased.This allows us to put quite stringent bounds on both these parameters.For 76 Ge nucleus, the lightest neutrino mass should be less than 10 −7 eV for θ LR ∼ 10 −4 while for a lightest neutrino mass of around 10 −3 eV the value of θ LR 10 −8 , where we consider a normal hierarchy among the active neutrino states.The bounds on 136 Xe nucleus are even more stringent with m ν1 10 −8 eV for θ LR ∼ 10 −4 and θ LR 10 −8 for m ν1 ∼ 10 −4 eV.Thus we can significantly constrain the model parameters in this case from the 0νββ studies.
For the Extended seesaw, we first consider one-generation scenario where in addition to one SM neutrino two sterile neutrinos are also present.Among the two sterile neutrinos one of them is very heavy with mass of 10 5 GeV leading to a negligible contribution in 0νββ process.The other sterile neutrino has a mass varying in between eV to MeV and this contributes significantly to the above mentioned process.We analyse a number of constraints on the model parameters, arising from 0νββ, reactor anti-neutrino experiment Daya Bay, non-unitarity constraint on the mixing matrix, as well as, theory-constraints.We further extend this simplistic one-generation analysis to higher generation with three active neutrinos and six sterile neutrinos for which we satisfy the neutrino oscillation data.We present a number of correlations between the mass of the lightest sterile neutrino, model parameters and several neutrino oscillation parameters.In threegeneration case, the mass of the heavy sterile states (N s) have been varied from 50 to 500 GeV.These sates give negligible contributions to 0νββ process, while the other three relatively light sterile states (S L s) give substantial contributions.With the considered parameter range we obtain the upper bound on the mass of the lightest sterile neutrino S 1 as 10 −6 GeV after imposing constraints from non-unitarity, 0νββ and others.The non-unitarity of PMNS matrix has direct impact on the mixing elements; in Extended seesaw, non-unitarity can be governed by the ratio of the bilinear mass term between active neutrino-light sterile states (M D ) to the corresponding bilinear mass terms among the sterile states (M S ).The upper bound on the ratio M D /M S is ∼ 10 −1 .We choose the effective neutrino mass scale (µM 2 D /M 2 S ) to be of the order of 10 −3 eV and we conclude in this scenario the lower bound on the elements of the complex symmetric matrix µ to be of the order of 10 −9 GeV.Another important constraint in our scenario is m s > µ below which the seesaw approximation ceases to be valid.It evidently gives lower bound on m si ; i.e., µ ii > 10 −9 GeV implies m si > 1 eV.The masses of the other two sterile neutrinos vary as m s2,3 ∼ O(1 − 10) eV.

Figure 1 :
Figure 1: Neutrino Majorana mass generation at one-loop in the LRS Zee model.

Figure 3 :
Figure 3: Feynman diagrams of all possible 0νββ processes in the LRS Zee model.

4 Figure 5 :
Figure 5: Half-life of 0νββ process for 76 Ge (left panel) and 136 Xe nucleus (right panel) as a function of left-right charged gauge boson mixing.

θLR = 10 Figure 6 :
Figure 6: Half-life of 0νββ process for 76 Ge and 136 Xe nucleus as a function of CP violating phase in the PMNS matrix.

Figure 7 :
Figure 7: Half-life of 0νββ process for 76 Ge nucleus as a function of the the lightest neutrino mass for several fixed values of θ LR and δ CP .

Figure 8 :
Figure 8: Half-life of 0νββ process for 76 Ge nucleus as a function of left-right charged gauge boson mixing for several fixed values of m ν1 and δ CP .

Figure 9 :m ν 1 igure 10 :
Figure 9: Half-life of 0νββ process for 76 Ge nucleus as a function of the CP violating phase of the PMNS matrix for several fixed values of θ LR and m ν1 .

Figure 11 :
Figure 11: Allowed parameter space of light sterile neutrino in (a) MeV and in (b) eV range as function of the model parameters M D and M S in Extended Seesaw scheme.The regions have been obtained from theoretical constraints, light neutrino measurements and 0νββ results.The parameter µ has been set to 10 −2 GeV and 10 −10 GeV for MeV and eV range respectively.In both cases, M R = 10 5GeV.The orange shaded region of (b) will be constrained further from reactor anti-neutrino Daya Bay experiment[97] that will be elaborated in Fig.12.
of M S − M D plane in ∆m 2 s1 (eV 2 ) − θ 21s in the Fig.12.The filled-in yellow box covered by the magenta color line corresponds to the aforementioned region of eV sterile of Fig.11b.The lower line corresponds to the value M S = 0.0032 GeV where as upper line corresponds to M S = 0.01 GeV.
, we have shown the scatter plot in plane M D22 −M S22 and µ I 23 (GeV) − M D11

23 10 − 9
compared to the rest of the region.
each nucleus ( 76 Ge or 136 Xe), we thus get four separate cases which are: (a) H min − ( 76 Ge/ 136 Xe) min : Corresponds to the case where the Higgs boson mixing is minimum and minimum value for the 76 Ge/ 136 Xe NME has been used.(b) H max − ( 76 Ge/ 136 Xe) min : Corresponds to the case where the Higgs boson mixing is maximum and minimum value for the 76 Ge/ 136 Xe NME has been used.(c) H min − ( 76 Ge/ 136 Xe) max : Corresponds to the case where the Higgs boson mixing is minimum and maximum value for the 76 Ge/ 136 Xe NME has been used.(d) H max − ( 76 Ge/ 136 Xe) max : Corresponds to the case where the Higgs boson mixing is maximum and maximum value for the 76 Ge/ 136 Xe NME has been used.