NMEs for $0\nu\beta\beta(0^{+}\rightarrow 2^{+})$ of two-nucleon mechanism for ⁷⁶Ge

In the standard model, the nuclear weak decay is interpreted as the low-energy effective theory for weak interaction. This decay is mediated by the left-handed gauge boson $W^{\pm}$. The mass of $W^{\pm}$ are acquired through the spontaneous symmetry breaking by the so-called Higgs mechanism. However, the Yukawa coupling of Higgs particle to neutrinos is absent in the standard model due to the absence of the right-handed neutrinos. The discovery of neutrino masses from oscillation experiments then asks for new physics beyond the standard model. As an extension to Standard model, the L-R symmetric model PS74 ; MP75 ; SM75 introduces the right-handed $SU(2)_{R}$ gauge symmetry and a hence heavy right-handed gauge boson from symmetry breaking with extra Higgs bosons at a higher energy scale beyond electroweak scale. In such a theory, the introduction of lepton number violating neutrino Majorana mass terms together with normal Dirac mass terms gives naturally the tiny neutrino mass through the so-called See-Saw mechanism Kin03 . Such extensions to the Standard Model could also affect the rare nuclear process called neutrinoless double beta decay ($0\nu\beta\beta$). The participation of right-handed weak gauge bosons will also induce the emission of right handed leptons. The simultaneous presence of weak currents with both chirality will introduce a momentum term into the neutrino propagator. These terms are not suppressed like the mass terms due to the smallness of neutrino mass. The right-handed weak currents will contribute to the decay to the ground states with extra terms and change the electron spectra SDS15 . Nevertheless, these terms are suppressed by the new physics parameters as well as the electron wave functions for p partial waves. Therefore, they are hindered in normal neutrinoless double beta decay compared to the neutrino mass mechanisms. On the other hand, the decay to the $2^{+}$ states, are dominated by the helicity changing mechanisms (V+A terms, see DKT85 ). In this sense, the branching ratio of neutrinoless double beta decay to the $2^{+}$ state (hereafter $0\nu\beta\beta(2^{+})$, the spin-parity of the final states of the decay are included inside the parenthesis ) could help to reveal the underlying mechanisms of this very rare decay. Nevertheless, experimentally such a process is extremely difficult to observe due to the large $2\nu\beta\beta(0^{+})$ background around the position of $0\nu\beta\beta(2^{+})$ Q value. Despite the difficulties, the observation of $0\nu\beta\beta(2^{+})$ together with that for the decay to the ground states will determine the underlying mechanisms of the neutrinoless double beta decay and pave our way to new physics beyond the standard model. For example, the observation of $0\nu\beta\beta(2^{+})$ could possibly rule out a category of mechanisms where no right-handed gauge bosons or fermions are present.

There are numerous publications dedicated to the nuclear many-body calculations for neutrinoless double beta decay with various approaches, e. g. the Shell Model calculations CNP99 ; Men17 , the QRPA calculations SPV99 ; SRF13 , the IBM calculations BKI13 and the nuclear meanfield calculations SYR14 , especially the recently developed ab initio methods YBE19 . However all these works focus on $0\nu\beta\beta(0^{+})$, there aren’t to many theoretical investigations available for $0\nu\beta\beta(2^{+})$ in the literature for the matrix elements (NME). An earlier calculation Tom88 with the Projected Hartree Fock Bogoliubov method (PHFB) suggests, that the nuclear matrix elements (NME) for this decay mode are much smaller than for the decay to ground states. The QRPA method is widely used in double beta decay calculations SRF13 ; ME13 ; HS15 . The QRPA can also be adopted to describe the vibrational $2^{+}$ states. Therefore, there are attempts to use QRPA to calculate $\beta\beta$-decays to the first $2^{+}$ state SC93 ; SSF97 . And most of them focus on the $2\nu\beta\beta$ case. In this study, we go step forward by carrying out the QRPA calculations for $0\nu\beta\beta(2^{+})$ for ⁷⁶Ge. Taking advantage of the QRPA method, we take into account the contributions from all the intermediate states. Also we can include the isoscalar particle-particle residual interaction which is missing in PHFB calculations. At this first attempt, we do not include too many examples. We focus on one nucleus, ⁷⁶Ge which is also the candidate treated in ref. Tom88 . It has been shown in Tom88 , that besides the nucleon mechanism, the N* could also play an important role, we will not discuss this in the current study. Also, as in CNP99 , the induced weak current will further reduce the NME, this will be neglected in this work. As well as the heavily suppressed neutrino mass mechanism for $0\nu\beta\beta(2^{+})$ through nuclear recoil Tom99 .

The current article is arranged as follows: first we give the general formalism of our many-body calculations, and then we show the results and discuss possible uncertainties, and finally we present the conclusions and outlook.

The half-lives of $0\nu\beta\beta(2^{+})$ can be expressed in a simple form, while we consider the light neutrino only DKT85 ; Tom88 :

Here $F_{1(2)}$ are the phase space factors expressed in DKT85 ; Tom88 . $\langle\lambda\rangle$ and $\langle\eta\rangle$ are the new physics parameters which are model dependent. In the L-R symmetric model MP75 : $SU(2)_{L}\times SU(2)_{R}\times U(1)_{B-L}$, we have DKT85 :

Here $U_{ej}$ and $V_{ej}$ are matrix elements for the generalized PMNS matrix Xin11 . $\lambda\approx(M_{WL}/M_{WR})^{2}$ are the square of the ratio of the masses between the mass eigenvalues of the light and heavy gauge bosons, $\eta\approx\tan\xi$ is the mixing angle between the left-handed gauge boson and the heavy gauge boson mass eigenvalues.

$M_{\lambda}$, $M_{\eta}$ and $M^{\prime}_{\eta}$ are the nuclear matrix elments (NMEs) as combinations of different components DKT85 :

The different coefficients $C_{Ii}$ are given in table I of Tom88 .

The general form of above NMEs can be expressed as:

A general derivation of the single particle matrix elements are given in the appendix. For the related operators, we have the form Tom88 :

And:

Here $r$ is the relative distance between the two decaying nucleons $\vec{r}=\vec{r}_{1}-\vec{r}_{2}$ and $\vec{r}_{+}=(\vec{r}_{1}+\vec{r}_{2})/2$ is their average position.

Here, $M_{1}$, $M_{2}$ and $M_{3}$ are the space-space components of the current-current interaction, while $M_{4}$ is the time-time components and $M_{5}$, $M_{6}$ and $M_{7}$ are the time-space components. These time-space components of the NME’s appear only in the L-R symmetric case and are missing in the neutrino mass mechanisms. In all these NME’s, we find a similarity between $M_{2}$ and $M_{GT}^{0\nu}$ for $0\nu\beta\beta(0^{+})$ as well as $M_{3}$ and $M_{T}^{0\nu}$. The GT operator $\sigma$ or the tensor operator $[\sigma\otimes\hat{r}]^{2}$ replace the scalar products in $0\nu\beta\beta(0^{+})$. We also find an analog similarity between $M_{4}$ and $M_{F}^{0\nu}$, where $\hat{r}$’s come from the p-wave electron and the momentum term form a tensor product instead of a scalar product in $0\nu\beta\beta(0^{+})$.

The neutrino potential differs from that of the mass mechanism due to the momentum terms in the neutrino propagator SDS15 :

By deriving this, we assume that the two electrons share the decay energy, therefore $E_{N}=E_{m}+M_{m}-(M_{i}+M_{f})/2$ and $E_{m}$ is the excitation energy of the m$th$ excited states of the intermediate nucleus. The nuclear radius $R=1.2A^{1/3}\b{[}fm]$ introduced here makes the final NME dimensionless. For the form factor $F(q^{2})$ we use a dipole form with the parameters are as in SPV99 .

Usually, an extra radial function $f(r)$ should be multiplied to the above expression to take into consideration the strong repulsive nature of nucleon-nucleon interaction at short range. This is usually called the short range correlation (src) function, and in our calculation we choose the CD-Bonn or Argonne src extracted from the corresponding nuclear force with the form given in SFM09 .

For the reduced one-body density for transitions from intermediate states to final $2^{+}$, the expression is complicated FF20 :

Here $X$’s and $Y$’s are the amplitudes for pn-QRPA (proton-neutron Quasi-particle Random Phase Approximation) describing the intermediate states and $\mathcal{X}$’s and $\mathcal{Y}$’s are the amplitudes for CC-QRPA (Charge Conserving QRPA) describing the final $2^{+}$ state FF20 . And $u$’s and $v$’s are the BCS coefficients.

The reduced one body density for transitions from the initial states to the intermediate states can be expressed as SPV99 :

We also introduce the overlap between the initial and final intermediate states with the form SPV99 :

For simplicity, we set ${}_{f}\langle BCS|BCS\rangle_{i}\approx 1$.

The details of derivations of the BCS coefficients and respective QRPA amplitudes for the current work are presented in FF20 .

For our QRPA calculations, the single particle energies are taken from the solutions of Schrödinger equations with a Coulomb corrected Woods-Saxon potential. For the single particle wave functions, we use the Harmonic Oscillator wave functions. For the pairing part, we use the realistic CD-Bonn force derived from Brückner G-matrix. This is also used for the pn-QRPA and CC-QRPA residual interactions. A fine tuning of the interactions is needed to reproduce the experimental values. For the pairing part, we fit the two parameters $g_{pair}^{p}$ and $g_{pair}^{n}$ to reproduce the odd-even mass staggering. For pn-QRPA, we multiply the G-matrix by overall renormalization factors $g_{ph}$ and $g_{pp}$’s for particle-hole and particle-particle parts, respectively. We set $g_{ph}=1$. And for $g_{pp}$, we fit the iso-scalar channel ($g_{pp}^{T=0}$) and iso-vector channel ($g_{pp}^{T=1}$) separately. $g_{pp}^{T=0}$ is fixed by reproducing the experimental $2\nu\beta\beta$ and $g_{pp}^{T=1}$ are fixed to put to zero the $2\nu\beta\beta$ Fermi matrix elements due to isospin symmetry restoration RF11 . For more details, one can refer to the our previous work FF20 . For a baseline calculation, we consider the CD-Bonn src bare axial vector coupling constant $g_{A0}=1.27$. And we use a model space with $N_{max}=6$ which consists of 28 single particle levels for both neutrons and protons.

In table.1, we present the NMEs for each operator. As a comparison, we also present the results from PHFB calculations. The current results (the baseline results) differ from the PHFB results by factors from five to more than one order of magnitude case by case. The largest deviation we see in $M_{7}$ for the $M^{\prime}_{\eta}$ part. We obtain $M_{7}$ much larger than $M_{6}$ in magnitude. Our results have have also different phases for these two NMEs, this then combined with the $C^{\prime}_{\eta}$ coefficients leads to the enhancement instead of the cancellations of $M^{\prime}_{\eta}$. Therefore, our $M^{\prime}_{\eta}$ is much larger than a previous PHFB calculation. PHFB gets an approximate cancellation between $M_{6}$ and $M_{7}$, which leads to an negligible $M^{\prime}_{\eta}$.

For the $M_{\lambda(\eta)}$ part, we find the NME’s have basically similar phases as the PHFB results except for $M_{2}$. On the other hand, our results are about one order of magnitude larger, although the relative ratios among different NMEs ($M_{1}-M_{5}$) are similar. Of these NMEs, $M_{1}$ is the largest. The second largest is $M_{2}$ and third is $M_{4}$. $M_{3}$ and $M_{5}$ are relatively small and hence less important. For $M_{\lambda}$, if we multiply the NMEs with the corresponding $C_{\lambda}$’s, we find that $M_{1}$ and $M_{2}$ contributes coherently, they are then cancelled by $M_{4}$, while the rest two NME’s contributes less than 10%. For $M_{\eta}$, all these three NMEs gives additive contributions, this makes $M_{\eta}$ about three times larger than $M_{\lambda}$. This is the reason, why $M_{\eta}$ is larger than $M_{\lambda}$ as also observed in Tom88 . In our case, this ratio is about three. This is similar to PHFB calculations, however, in their calculation, the strong cancellation gives a negligible $M_{\lambda}$, the ratio is about 30. In short, in their calculation, only $M_{\eta}$ is important while $M_{\lambda}$ and $M^{\prime}_{\eta}$ can be neglected due to the cancellations between different parts of the NME. We get quite different results, and $M^{\prime}_{\eta}$ is the most important contribution, with a value about three times larger than $M_{\eta}$. This will affect the constraints on new physics models and needs further investigation.

As we show above, no obvious suppression of $0\nu\beta\beta(2^{+})$ NME’s as claimed by Tom88 is found from current calculations compared to $0\nu\beta\beta(0^{+})$, this also agrees with the $q$ terms in $0\nu\beta\beta(0^{+})$ calculations MBK89 . This is the major difference of the current work and Tom88 . This is also different from $2\nu\beta\beta(2^{+})$, where the NME is suppressed by the cubic dependence of the energy denominator DKT85 . In this sense, suppression of $0\nu\beta\beta(2^{+})$, if it exists, must be related to other issues. This may come from the uncertainties of the many-body approaches, such as the size of the model space or other structure ingredients for the $2^{+}$ states which may lead to different transition rates from the intermediate states for various transitions. A more thorough comparative study could give us mor detailed hints.

Compared to previous calculations with PHFB Tom88 , the QRPA calculation goes beyond the closure approximation. We calculate explicitly the contribution from each intermediate state. In fig. (1-4), we present the individual contributions from different multipoles of the intermediate states and we will also show how the different approximations may affect these results. In each graph, the results are compared with our present baseline calculations with standard conditions described above.

For details of the structure of the nuclear MME’s, we start with our baseline calculation (e.g. orange bars in fig.1). For $M_{1}$, the multipoles give positive contributions with several exceptions. Unlike $0\nu\beta\beta(0^{+})$ ($M_{GT}^{0\nu}$), where the largest contributions come from low-spin intermediate states and the NME values decrease as spins increase, $M_{1}$ has its largest contribution from $4^{-}$. We find a rough trend that the NME’s first increase and then decrease as spin increases. And as one would expect, the NME’s from states with very high spin can be safely neglected. $M_{2}$ has basically the same characters as $M_{1}$ except the much smaller magnitude. A large contribution from $1^{-}$ is observed for $M_{2}$ but not for $M_{1}$. For most multipoles, $M_{2}$ has different signs as $M_{1}$, this contradicts conclusions in Tom88 . Not all multipoles contribute equally with similar spins, we find that the states with negative parity generally contribute more. In some sense, these two NME’s behave like $M^{0\nu}_{GT}$ for $0\nu\beta\beta(0^{+})$ as we mentioned above. The smallness of final $M_{3}$ comes partially from its magnitude and partially from the cancellations between low-spin and high-spin multipoles. This is analog to $M^{0\nu}_{T}$, although these two NME’s depend differently on $\hat{r}$.

$M_{4}$, the time-time component of the NME is on the other hand, very close to $M_{F}^{0\nu}$. The two NME’s have one thing in common: Only intermediate states with natural parity ($\pi=(-1)^{J}$) have non-zero contributions. In the current calculation, all multipoles contribute negatively except $1^{-}$. All these multipoles contribute basically with the same magnitude as $M_{2}$ and $M_{4}$ is the major cancellation to $M_{\lambda}$.

For the space-time components, $M_{5}$, $M_{6}$ and $M_{7}$, we find a quite different behavior. $M_{5}$ from each multipole is about one order of magnitude smaller than $M_{6}$ and $M_{7}$. Due to the strong cancellations among different multipoles, $M_{5}$ is one of the smallest of all the NME’s. Almost all multipoles contribute positively except a strong cancellation from $3^{-}$. For $M_{7}$, the important contributions comes from three multipoles ($3^{-}$, $5^{-}$ and $7^{-}$), all the other multipoles contribute much less. The lack of cancellations from other multipoles thus makes $M_{7}$ the largest from all the NME’s.

A possible cause of the smallness of NMEs in ref. Tom88 may come from the small model space used whith only two major shells. To test this, we plot in fig.1 the results with $N_{max}=5$ (blue bars) and $N_{max}=7$ (green bars). The sensitivity of the NME’s to diffeent model spaces are different. Also different is the sensitivity of the individual multipoles for each NME. For $M_{1}$, the influence from model space is generally small, there is no unique trend for different multipoles under the change of the model space. For some multipoles, the NME decreases with enlarged model space, but most cases we find that the extra orbitals will first enhance but then reduce the NME. The orbitals of different parity contribute to the NME differently. The addition of N=6 shell brings in the positive parity orbitals which enhance the NMEs for multipoles such as $3^{-}$ or $5^{+}$, while the negative parity orbitals from N=7 shells reduce the NME’s. Unlike the case of $M_{1}$, the increase of the model space enlarge $M_{2}$, especially for $1^{+}$, $2^{-}$ etc. For $1^{+}$, a strong suppression is observed when the N=6 shell is added. But such addition gives strong enhancement to $2^{-}$ or $3^{+}$. The addition of the N=7 shell to the model space causes milder changes. We see enhancement from $1^{-}$ and $2^{-}$ but reductions from $3^{-}$ and $4^{-}$ intermediate states. A similar behavior shows the $M_{3}$, where we find a strong enhancement from $1^{+}$ too. For all other multipoles, the change due to model space enlargement is relatively small. For $M_{4}$, low spin multipoles are much more sensitive to the model space changes and different multipoles behave differently, though we cannot find any specific patterns. This is also true for the space-time components of the NME. In general, they are less affected by the change of model space in our calculations.

If we look at the total changes of each NME from $N_{max}=5$ to $N_{max}=6$, $M_{1}$ increases about $10\%$. this is the smallest change among all the NMEs. Meanwhile, all the other NMEs changes about $20\%$ or more. By percentage $M_{3}$ changes by more than $200\%$. By the absolute values of NMEs, $M_{1}$, $M_{2}$, $M_{3}$ and $M_{4}$ get enhanced while the rest get reduced. When the $N=7$ shell is added, the change is relatively milder, especially for $M_{2}$, $M_{3}$ and $M_{4}$, this implies a general trend of convergence of the results with a larger model space. But for $M_{1}$, $M_{5}$, $M_{6}$ and $M_{7}$, we find a slower trend of convergence than for the above NMEs. In either cases, the changes from adding the $N=7$ shell is smaller than adding the $N=6$ shell. Current results suggest that the errors of adopting the current model space are generally smaller than $20\%$. Also, these results suggest that the deviation of our calculations and those in Tom88 is not caused by the small model space, which they adopted. It is most probable that the smallness of their results are caused by different nuclear structures. Additional work is needed to clarify this.

The form factor with a dipole form is widely used in $0\nu\beta\beta(0^{+})$ calculations SPV99 . In our calculations, we use the same form for $g_{V}(q^{2})$ and $g_{A}(q^{2})$. We find that the form factors are not very important except for the $1^{+}$ states of $M_{2}$. This may suggest that with the actual neutrino potential the low momentum parts where the $q^{2}$ satisfies $g_{A}(q^{2})\sim g_{A}(0)$ dominate, and the high momentum parts are either small or cancels each other. A careful check suggests the latter should apply. These behaviors also helps to explain the large reduction for some NME’s when the realistic neutrino potential is considered. With low momenta, the intermediate energies become important. In this sense, the choice of the intermediate energies becomes important. For calculations with closure approximation, a closure energy is needed. So we also check how the use of closure energy would affect the NME’s. This is illustrated with the red bars in fig.2. A closure energy of $7$ MeV is used. We find for most multipoles of most NME’s, the proper choice of closure energy brings small errors to the calculation and we can draw the conclusion, that using the closure energy in $0\nu\beta\beta(2^{+})$ barely changes the final results, the errors are within several percents, this agrees with $0\nu\beta\beta(0^{+})$ QRPA calculations. Another important correction comes from the induced weak current SPV99 , it is not included in the current calculations and will be implemented in our future study. If that is included, the NME will most probably be further reduced and new potentials needs to be introduced SDS15 .

In $0\nu\beta\beta(0^{+})$, src only plays an important role for the heavy neutrino mass mechanism FFS18 ; SYR14 , while for the light neutrino mass mechanism, the correction is relatively small, up to only several percents FFS18 ; SYR14 . In current calculations, we adopt two src’s SFM09 : the CD-Bonn type and Argonne type, they are obtained by fitting the respective nuclear potentials. From fig.3, we find similar trends as $0\nu\beta\beta(0^{+})$, the general correction from src is about several percent and for most cases, the two src’s give similar amount of corrections. For almost all multipoles, we find that the src enhances the NMEs more or less. However, the net effects to the NMEs are slight reductions for $M_{3}$ and $M_{5}$ due to cancellations among different multipoles. And enhancements of NME values for other operators are observed. The largest correction comes from $M_{6}$, then $M_{7}$ and $M_{2}$. For $M_{6}$ we can also find slight difference between the two src’s.

The most important parameters in QRPA calculations are the particle-particle interaction strength $g_{pp}$’s. For $0\nu\beta\beta(0^{+})$, one finds that $M_{GT}$ for both $2\nu\beta\beta$ and $0\nu\beta\beta$ depends sensitively on the isoscalar strength $g_{pp}^{T=0}$ while $M_{F}^{0\nu}$ are sensitive to the isovector strength $g_{pp}^{T=1}$. So we also test such dependence in current calculations. As we have shown in FF20 for the case of $2\nu\beta\beta$ to $2^{+}$ states, the GT type decays are also sensitive to $g_{pp}^{T=1}$ since the isovector particle-particle residue force affects the structure of $2^{+}$ states in QRPA calculations. In fig.4, we first switch off the isoscalar interaction (the blue bar). Then the effects of this interaction to NMEs can be estimated by comparing the blue and green bars. The results suggest that, while $M_{6}$ and $M_{7}$ are not closely related to $g_{pp}^{T=0}$, $M_{2}$ comes out to be the most sensitive one like its $0\nu\beta\beta(0^{+})$ counterpart $M_{GT}^{0\nu}$. Similar as $M_{GT}^{0\nu}$, $1^{+}$ comes out to be the most sensitive multipole, the increase of the $g_{pp}^{T=0}$ drastically changes the NME most probably due to the SU(4) symmetry restorationRF11 . And the effects of isoscalar residue interactions to the specific multipoles of specific NMEs are quite different, some NMEs are enhanced and some are reduced. In total, the introduction of isoscalar interactions enhances $M_{1}$, $M_{2}$ and $M_{7}$, but reduces other NMEs. And the magnitudes of these enhancements and reductions are really case dependent. The similar thing happens to isovector interactions (this interaction has been switched off for the red bars in fig.4), for some cases they are more important than the isoscalar interactions. In general, the particle-particle interaction is one of the most important source of the errors for QRPA calculations. $M_{4}$ is at least sensitive to $g_{pp}$’s while $M_{2}$ and $M_{7}$ are really sensitive to $g_{pp}^{T=0}$ and $g_{pp}^{T=1}$ respectively. Especially $M_{2}$, the absence of isoscalar interactions will reduce the NME by more than $30\%$.

The special attention should be paid to the case of the quenched $g_{A}$. In all our above analysis, we assume that $g_{A}$ is not quenched, however, in nuclear medium quenching of $g_{A}$ is observed. In current calculations, we use the simplest treatment for quenching, that is simply change the value of $g_{A}$, while there exists more fundamental approaches for the quenching of $g_{A}$ using the chiral two-body currentsESV14 . In our case, the difference between quenched $g_{A}$ calculations and our baseline calculations are the different fit of $g_{pp}^{T=0}$. Therefore, the difference for the individual NMEs are small as in fig.4. But the quenched $g_{A}$ will also change the coefficients $C$’s. As a result, the NMEs $M_{\lambda}$, $M_{\eta}$ and $M^{\prime}_{\eta}$ have been changed too. In general, these NME’s are reduced due to the convention used inTom88 . For $M^{\prime}_{\eta}$ a reduction of a rough factor of $g_{A}/g_{A0}$ is expected, since the two components $M_{6}$ and $M_{7}$ has the same dependence on $g_{A}$. For $M_{\lambda}$ and $M_{\eta}$, different components have a different $g_{A}$ dependence. And if we take $g_{A}=0.75g_{A0}$ we find a rough cancellation for $M_{\lambda}$ as predicted in Tom88 , this comes from the interplay among different components. Meanwhile the reduction for $M_{\eta}$ is basically $25\%$. The severe reductions of $M_{\lambda}$ emphasize the importance of where the quenching of $g_{A}$ originates and how to treat it properly.

In this study, we calculate the NME of $0\nu\beta\beta$ to $2^{+}_{1}$ for ⁷⁶Ge. We got quite large results compared to previous calculations. We estimate the errors of current calculation by changing several parameters we use. We find that $g_{A}$ may be a very important issue for the final NME’s. Further investigations, such as the role of the induced weak current, the anharmonicity beyond QRPA, and the decay mechanism mediated by N*, are needed for much more detailed conclusions.

This work is supported by the ”Light of West China” program and the ”From 0 to 1 innovative research” program both from CAS.