Interacting Shell Model Calculations for Neutrinoless Double Beta Decay of 82Se With Left-Right Weak Boson Exchange

In the present work, the λ mechanism (left-right weak boson exchange) and the light neutrino-exchange mechanism of neutrinoless double beta decay is studied. In particular, much attention is paid to the calculation of nuclear matrix elements for one of the neutrinoless double beta decaying isotopes 82Se. The interacting shell model framework is used to calculate the nuclear matrix element. The widely used closure approximation is adopted. The higher-order effect of the pseudoscalar term of nucleon current is also included in some of the nuclear matrix elements that result in larger Gamow-Teller matrix elements for the λ mechanism. Bounds on Majorana neutrino mass and lepton number violating parameters are also derived using the calculated nuclear matrix elements.


INTRODUCTION
Neutrinoless double beta decay (0]ββ) is a rare second-order weak nuclear process. In this process, neutrino comes as a virtual intermediate particle when two neutron pairs decay into two proton pairs inside some even-even nuclei. Thus, it violates the lepton number by two units. The 0]ββ experiment is one of the possible ways to determine the effective neutrino mass (Schechter and Valle, 1982;Tomoda, 1991;Avignone et al., 2008;Rodejohann, 2011;Deppisch et al., 2012) and can help to solve many mysteries of neutrinos, such as whether neutrinos are their own anti-particle (Majorana neutrino) or not (Dirac neutrino) (Schechter and Valle, 1982;Rodejohann, 2011;Deppisch et al., 2012).
In the present work, we focus on the left-right weak boson (W L -W R ) exchange λ mechanism along with the standard light neutrinoexchange mechanism (W L − W L exchange) of the 0]ββ mediated by light neutrinos (Bhupal Dev et al., 2015;Horoi and Neacsu, 2016;Šimkovic et al., 2017). The λ mechanism has origin in the left-right symmetric mechanism with right-handed gauge boson at the TeV scale (Šimkovic et al., 2017). Thus, it will be interesting to study how the λ mechanism can compete with the standard light neutrinoexchange mechanism when both the mechanisms co-exist. Hence, in the present work, we are eager to study the λ and light neutrinoexchange mechanisms together.
In left-right symmetric model, there is another mass independent mechanism called η mechanism which occurs through W L − W R mixing. It will be interesting to study η mechanism along with λ mechanism of 0]ββ. But, η mechanism is suppressed due to W L − W R mixing as compared to λ mechanism (Barry and Rodejohann, 2013). Hence, in the present work, we are interested to study the mass independent λ mechanism along with the mass dependent standard light neutrino-exchange mechanism. In future studies, we will extensively explore the η mechanism of 0]ββ along with other mass independent and dependent mechanisms in left-right symmetric model.
One of the motivations of the present work is to include effects of some of the revisited formalism of Ref. (Štefánik et al., 2015) on light neutrino-exchange and λ mechanism of 0]ββ. The revised formalism was exploited to include the effects of the pseudoscalar term of nucleon currents. Using the revised formalism of Ref. (Štefánik et al., 2015), the NMEs for λ, and light neutrino-exchange mechanisms of 0]ββ are calculated using the QRPA model for several 0]ββ decaying isotopes using closure approximation in Ref. (Šimkovic et al., 2017). Most of the NMEs relevant for λ and light neutrino-exchange mechanisms are also calculated using ISM in Ref. (Horoi and Neacsu, 2018) using the closure approximation for different 0]ββ decaying isotopes (including 82 Se). In this case, some of the NMEs are calculated without including the higher-order terms (for example, pseudoscalar and weak magnetism terms) of the nucleon currents. Recently, using the revised formalism of Ref. (Štefánik et al., 2015), we have also calculated the NMEs for 48 Ca in Ref. (Sarkar et al., 2020a) using the non-closure approximation and found a significant change in some of the NMEs for including the pseudoscalar term. Thus, we have tried here to include the revised higher-order effect of the pseudoscalar term of nucleon current for the λ mechanism of 0]ββ of 82 Se using ISM. The 0]ββ of 82 Se is one of the experimental interests of CUPID (Dolinski et al., 2019;Pagnanini et al., 2019) and NEMO-3 (Arnold et al., 2020) experiments. Hence, it is important to study the nuclear structure aspects of 0]ββ of 82 Se theoretically. In recent years, one of the most important studies on light neutrinoexchange 0]ββ of 82 Se was performed in the ISM framework in Ref.  using the non-closure approximation. Here we focus on the λ mechanism of 0]ββ of 82 Se in the closure approximation using the revised nucleon current term.
This paper is organized as follows. In Section 2, the expression for decay rate and the theoretical formalism to calculate NMEs for the λ and light neutrino-exchange mechanisms of 0]ββ are presented. The results and discussion are presented in Section 3. A summary of the work is given in Section 4.

Decay Rate for λ Mechanism of 0νββ
If both light neutrino-exchange (W L − W L exchange) and λ mechanisms (W L − W R exchange) of 0]ββ co-exist, one can write the decay rate for 0]ββ as (Štefánik et al., 2015; Šimkovic et al., 2017) where the coupling constant λ is defined as (Šimkovic et al., 2017) λ The M W L and M W R are masses of the Standard Model lefthanded W L and right-handed W R gauge bosons, respectively. The η ] of Eq. 1 is an effective lepton number violating parameters for W L − W L exchange, η λ is an effective lepton number violating parameters for W L − W R exchange, and ψ denotes the CP violating phase. These parameters are given in Ref. (Šimkovic et al., 2017) as Here, m ββ is the effective Majorana neutrino mass defined by the neutrino mass eigenvalues m j and the neutrino mixing matrix elements U ej (Horoi and Stoica, 2010): The U, and T are the 3 × 3 block matrices in flavor space, which constitute a generalization of the Pontecorvo-Maki-Nakagawa-Sakata matrix, namely the 6 × 6 unitary neutrino mixing matrix (Štefánik et al., 2015; Šimkovic et al., 2017).
The amplitude of λ mechanism is given by (Bhupal Dev et al., 2015) where λ is defined earlier, G F is the Fermi constant for weak interaction, and q is the virtual Majorana neutrino momentum. The coefficients C I (I mm, mλ and λλ) of Eq. 1 are linear combinations of products of nuclear matrix elements and phasespace factors (Šimkovic et al., 2017).

Nuclear Matrix Elements for λ Mechanism of 0νββ
Matrix elements required in the expression of C I are (Šimkovic et al., 2017).
There are two approximations for calculating the NME, one is non-closure approximation and another is the widely used closure approximation. In non-closure approximation, the radial neutrino potential H α (r, E k ) has explicit dependence on energy of the intermediate state |k〉. In non-closure approximation, the radial neutrino potential for λ mechanism of 0]ββ are is given as integral over Majorana neutrino momentum q (Sen'kov and Horoi, 2013): where R is the radius of the parent nucleus, and the f α (q, r) factor (Appendix B) contains the form factors that incorporates the effects of finite nucleon size (FNS), and higher-order currents (HOC) of nucleons (Šimkovic et al., 1999), which is given in Appendix B of the manuscript. The E i and E f are the g.s. energy of the initial and final nucleus of the 0]ββ decay, respectively. The non-closure approximation is computationally very challenging, because in this approximation, the NME has explicit dependence on the energy of large numbers of virtual intermediate state |k〉 and calculating these states requires enormous computational power. Particularly, for higher mass region isotopes, some of the calculations are still beyond the reach of current generation's high-performance computers. Fortunately, the most of the contributions on NME of 0]ββ come from low lying energy states up to 10-12 MeV of the intermediate nucleus (Sen'kov and Horoi, 2013;Sarkar et al., 2020a). Thus, one can replace the effects of E k with a suitable constant energy called closure energy 〈E〉 without affecting the value of NME too much, and this approximation is known as closure approximation. In this approximation, one assumes (Sen'kov and Horoi, 2013) and the radial neutrino potential operator of Eq. 16 becomes In closure approximation, the 0]ββ decay operators defined in Eq. 15 become O GT,ωGT,qGT 12 The closure approximation is widely used in literature as it eliminates the complexity of calculating a large number of virtual intermediate states (Horoi and Stoica, 2010;Sen'kov and Horoi, 2013;Sarkar et al., 2020a). One can find suitable values of 〈E〉 using the method described in Ref. (Sarkar et al., 2020a), such that using closure approximation, one can get NME near to the nonclosure approximation.
In the calculation of the NME of 0]ββ, it is also necessary to take into account the effects of short-range correlations (SRC). A standard method to include SRC is via a phenomenological Jastrow-like function (Vogel, 2012;Menéndez et al., 2009;Šimkovic et al., 2009). Including SRC effect in the Jastrow approach, one can write the NME of 0]ββ defined in Eq. 14 as (Vogel, 2012) where Jastrow-type SRC function is defined as In literature, three different SRC prametrizationparameterization are used: Miller-Spencer, Charge-Dependent Bonn (CD-Bonn), and Argonne V18 (AV18) to parametrize a, b, and c (Horoi and Stoica, 2010). These parameters are chosen in such a way that the two-body wave function of two-body matrix elements (TBME) for 0]ββ are still normalized. The parameters a, b, and c in different SRC parametrizations are given in Table 1.
This approach of using a Jastrow-like function to include the effects of SRC is extensively used in Refs. (Menéndez et al., 2009;Horoi and Stoica, 2010;Neacsu et al., 2012).

The Closure Method of Nuclear Matrix Elements Calculation for 0νββ in ISM
The (M GT,ωGT,qGT ) (M F,ωF,qF ), and (M T,ωT,qT ) matrix elements of the scalar two-body transition operator O α 12 of 0]ββ can be expressed as the sum over the product of the two-body transition density (TBTD) and anti-symmetric two-body matrix elements (〈k 1 ′ , k 2 ′ , JT|O α 12 |k 1 , k 2 , JT〉 A ) : where, α (F, GT, T, ωF, ωGT, ωT, qF, qGT, qT), J is the coupled spin of two decaying neutrons or two final created protons, τ − is the isospin annihilation operator, A denotes that the two-body matrix elements (TBME) (Appendix A) are obtained using anti-symmetric two-nucleon wavefunctions, and k 1 stands for the set of spherical quantum numbers (n 1 ; l 1 ; j 1 ) (similar definition for k 2 , k 1 ′, k 2 ′). The |i〉 is 0 + ground state (g.s.) of the parent nucleus, and |f〉 is the 0 + g.s of the granddaughter nucleus.
The TBTD can be expressed as ) where, andÃ are the two particle creation and annihilation operator of rank J, respectively. Most of the available public shell model code does not provide the option to calculate TBTD directly. One of the ways is to calculate TBTD in terms of a large number of two nucleon transfer amplitudes (TNA), assuming 0]ββ decay occurs through (n − 2) channel . In (n − 2) channel of 0]ββ, the TNA are calculated with a large set of intermediate states |m〉 of the (n − 2) nucleons system, where n is the number of nucleons for the parent nucleus. In this approach, the TBTD in terms of TNA is expressed as ) where, TNA are given by Here, J m is the spin of the allowed states |m〉 of intermediate nuclei. J 0 is spin of |i〉 and |f〉. J m J when J 0 0 .

RESULTS AND DISCUSSION
We have used JUN45 effective shell model Hamiltonian (Honma et al., 2009) of fpg model space to calculate the relevant initial, intermediate, and final nuclear states for 0]ββ of 82 Se. In the fpg model space, valence nucleons can occupy the orbitals 0f 5/2 , 1p 3/2 , 1p 1/2 , and 0g 9/2 . For the 0]ββ decay of 82 Se through (n − 2) channel, the states of allowed spin-parity of 80 Se acts as intermediate states for TNA calculations. The nuclear shell model code KSHELL (Shimizu et al., 2019) was used in the calculation. For comparing some of the TNA values, NushellX@MSU (Brown and Rae, 2014) code was also used. In the present calculation, we have included the first 100 states of different allowed spin-parity of 80 Se in calculating the TNA. Earlier, it was found that considering around the first 50 states is enough to get the saturated value of NME, as the most dominating contributions come from the first few initial states Sarkar et al., 2020b).
We have adopted the widely used closure approximation with the closure energy 〈E〉 0.5 MeV. Earlier studies of Refs. (Sarkar et al., 2020a;Sarkar et al., 2020b) suggested that 〈E〉 0.5 MeV is a suitable value that is close to optimal closure energy and, thus, gives NME near to the NME in the non-closure approximation. The non-closure method can give the exact value of NME, but the present study is beyond the scope of studying it. But, according to earlier results (Sarkar et al., 2020a;Sarkar et al., 2020b), with 〈E〉 0.5 MeV, one can get NME in the closure approximation close to the NME in non-closure approximation (within 1% difference).
Different types of NMEs for light neutrino-exchange and λ mechanism of 0]ββ for 82 Se is shown in Table 2. Here, NMEs are calculated in different SRC parameterization schemes. All standard effects of FNS + HOC are taken care of in all calculations. It is found that the Gamow-Teller matrix elements dominate over Fermi and tensor type matrix elements. Also, it is found that the M qGT type matrix element associated with the λ mechanism is relatively large as compared to standard light neutrino-exchange Gamow-Teller matrix element M GT . This leads to the large value of total NME M 1+ for λ mechanism as compared to total NME M ] for light neutrino-exchange mechanism.
This increment of M qGT type of NME, which is obtained through the new revised expression of the nucleon currents of Ref. (Šimkovic et al., 2017), is surprisingly high. It is coming through the new revised expression of the nucleon currents of Ref. (Šimkovic et al., 2017) which includes the higher-order term (pseudoscalar) of the nucleon currents. In our calculation, Eq. 39 is used to calculate M qGT type NME using the revised formalism of nucleon currents of Refs. (Štefánik et al., 2015; Šimkovic et al., 2017).
An old equivalent expression of Eq. 39 is also found in Ref. (Horoi and Neacsu, 2018), which one can write using Eq. (A2c) and Eq. (A4b) of Ref. (Horoi and Neacsu, 2018) as Using this old value of f qGT (q, r), the M qGT type NME will be significantly smaller, as reported earlier.
Here we include the higher-order current effect of pseudoscalar term in Eq. 39 as suggested in Ref. (Šimkovic et al., 2017) which is enhancing the M qGT type NME as compared to standard M GT type NME. A similar type of enhancement in M qGT type NME was also found in our earlier study for 48 Ca (Sarkar et al., 2020a).
We have also decomposed the NME in terms of coupled spinparity (J π ) of two decaying neutrons and two created protons in the decay. Decomposed NME gives us a picture of the role of individual spin-parity on NME. The contribution of NMEs through different J π is shown in Figures 1-3 for different types of NME. Figure 1 examines the decomposition for M F,GT,T type matrix elements associated with light neutrinoexchange mechanism, where Figures 2, 3 examine the NME as function of J π for M ωF,ωGT,ωT and M qF,qGT,qT type NMEs, respectively, for λ and interference mechanism. All results are presented for AV18 SRC parameterization.   For all types of NMEs, the most dominating contribution comes from 0 + states and 2 + states. The pairing effect is in play for dominating even-J π contributions . The NME from 0 + and 2 + states has opposite signs and, thus, cancel the effects of each other. Other non-negligible contributions come through 4 + , 3 − , 5 − , and 7 − states. Now we will discuss how the calculated NMEs will help to determine the bounds on Majorana neutrino mass and various lepton number violating parameters, using the lower limit on the experimental half-life of the decay. The inverse of half-life for 0]ββ is given in Eq. 1. It is found that the half-life is influenced by the term C I (I mm, mλ, λλ), lepton number violating term η ] and η λ , which are unknown, and CP-violating phase ψ. The C I are defined in Eq 7 and (8), (9), which contains mainly phase space factors and relevant NMEs. To calculate C I , we have used the improved values of phase space factors calculated in Ref. (Štefánik et al., 2015), and for the NMEs, we have used the results of Table 2 using ISM.
The results for C I of light neutrino-exchange and λ mechanisms of 0]ββ decay of 82 Se and 48 Ca are presented in Table 3. Here, the results for 48 Ca are taken from our earlier work using the closure approximation on the λ mechanism (Sarkar et al., 2020a). It is found that values of C mm (light neutrinoexchange) and C λλ (λ mechanism) are similar in values, which shows the dominance of each of these mechanisms on 0]ββ halflife. The interference term (C mλ ) of both the mechanisms are relatively smaller, which shows the less importance of the interference mechanism.
We have also calculated the upper bound on unknown Majorana neutrino mass (m ββ ) and lepton number violating parameter: the right-handed current coupling strength η λ , using the experimental constraint on T 0]−exp 1/2 of Ref. (Arnold et al., 2005) for 82 Se and of Ref. (Arnold et al., 2016) for 48 Ca. The upper limits on m ββ and η λ are also presented in Table 3 for 82 Se and 48 Ca when both light neutrino-exchange and λ mechanisms co-exist. With the experimental lower limit on T 0]−exp 1/2 , the upper limits on Majorana neutrino mass (m ββ ) are found to be 1.83 and 17.92 eV, respectively, for 82 Se and 48 Ca. This difference of m ββ value for 82 Se and 48 Ca is quite large and also found in earlier work (Šimkovic et al., 2017). With the recent progress and future prospects of new generation experiments, lower limits on T 0]−exp 1/2 will be gradually improved and thus, will improve the upper limit on m ββ and also reduce the differences for different isotopes.

SUMMARY
In summary, we have studied how the left-right weak boson exchange (λ) mechanism of 0]ββ decay is competing with the standard light neutrino-exchange mechanism. Our interest of isotope was one of the prominent 0]ββ decaying isotope 82 Se. Particularly, we have calculated the NMEs for 0]ββ of 82 Se when both standard light neutrino-exchange and λ mechanisms coexist. The revised formalism for nucleon currents to include the pseudoscalar term was taken care of. The nuclear shell model framework was used in the calculation, and the widely used closure approximation was adopted with suitable closure energy. Nuclear states of initial, final, and intermediate states are calculated for fpg model space with JUN45 effective shell model Hamiltonian using shell model code KSHELL. These nuclear states are used to calculate TNA, which comes in the expression of NME of 0]ββ through (n − 2) decay channel. Using the calculated NMEs, we have also calculated the upper bounds on Majorana neutrino mass and lepton number violating parameters.  is taken from the experimental lower limit on half-life from Ref. (Arnold et al., 2005) for 82 Se and from Ref. (Arnold et al., 2016) for 48 Ca. All results are for AV18 type SRC parameterizaionparameterization. We have assumed CP conservation (ψ 0). The results are compared with QRPA calculations for λ mechanism of Ref. (Šimkovic et al., 2017). The results show that particularly M qGT type matrix element of λ mechanism is significantly enhanced as compared to standard M GT type NME for the inclusion of the higher-order effect of the pseudoscalar term in the nucleon current. A similar type of enhancement in M qGT type NME was also found in our earlier study for 48 Ca (Sarkar et al., 2020a). The dominance of 0 + and 2 + states of neutron-neutron (proton-proton) pairs were also observed, just like earlier studies.
With the experimental lower limits on the half-life, we have used our calculated NMEs to set the upper bounds on Majorana neutrino mass (m ββ ). The upper limits of values of m ββ are found to be 1.83 and 17.92 eV, respectively, for 82 Se and 48 Ca. With the new generation of experiments, the lower limit on half-life will be further improved, and thus we can expect a much more refined upper bound on m ββ , which may be below 1 eV. Also, the difference for the value of m ββ will be reduced.
The term C I (I mm, mλ, λλ), which contains the phase space factors and NMEs, was also evaluated. The C mm for light neutrino exchange and C λλ for λ mechanism were found to be similar in values, that were larger than the term C mλ for the interference of both the mechanisms. This shows the dominance of light neutrino exchange and the λ mechanisms over the interference mechanism. The overall dominant effect of light neutrinoexchange mechanism is observed over λ mechanism and interference of both the mechanisms for very small values of lepton number violating η λ parameter.
In the future, it will be interesting to see the competing effect of the λ mechanism on the light neutrino-exchange mechanism and also how their contribution on 0]ββ half-life will be evaluated in the current and future generation experiments.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
The idea of the article was originated by YI. He has also contributed to calculating the nuclear states, interpretation of results, and manuscript writing. SS is responsible for the calculation of the TNA and TBME part of the NME. He has also actively participated in preparing the manuscript. Overall, both the authors have contributed enough in to various stages of preparing the final manuscript.

FUNDING
YI is grateful for the funding support from JSPS KAKENHI Grant No.17K05440.

ACKNOWLEDGMENTS
Numerical computation in this work was carried out at the Yukawa Institute Computer Facility. YI acknowledges the Tokyo Institute of Technology for allowing to use of the highperformance computing facility to perform nuclear states calculation using KSHELL. YI is also grateful to Prof. Noritaka Shimizu, CNS, the University of Tokyo, for providing the 2020 version of shell model code KSHELL.