Nuclear structure and double beta decay

Study of the neutrinoless double beta decay is the most sensitive current test of the total lepton number conservation. If observed, it would serve as a proof that not only the lepton number L_tot is not conserved, but also that neutrinos are massive Majorana fermions. Hence the study of this process is one of the important parts of the search for 'physics beyond the standard model'. Different aspects of the complex of problems associated with the ββ decay are described in this focus issue. Since the process involves nuclei, i.e. complicated many body systems, the analysis of the decay necessary involves various nuclear structure problems. In this paper the general questions of nuclear structure, relevant for the understanding of the ββ decay rate, are discussed. In other works in this issue the particular approximate techniques for the treatment of the corresponding nuclear matrix elements are described. Here I try to discuss the more general framework. Various aspects, both theoretical and experimental, of the ββ decay have been reviewed many times. I quote here just some of the review articles [1–5], earlier references can be found there.

In double beta decay two neutrons, bound in the ground state of an even–even initial (or parent) nucleus are transformed into two bound protons, typically also in the ground state of the final (or granddaughter) even–even nucleus, with the simultaneous emission of two electrons only for the 0νββ mode, or two electrons plus two $\bar{\nu }_e$ for the 2νββ mode. Transitions leading to the excited bound states of the final nucleus are sometimes kinematically allowed as well; however, we will concentrate here on the most often studied case of the ground state to ground state decays

$$\begin{equation} (Z,A)_{{\rm g.s.}} \rightarrow (Z+2,A)_{{\rm g.s.}} + 2e^- + (2\bar{\nu }_e) . \end{equation} \tag{ 1 }$$

Note that transitions of two bound protons into two bound neutrons (nuclear charge is decreased by two units in that case) with the emission of either two positrons, or a positron accompanied by the electron capture, or by two simultaneous electron captures, are also kinematically allowed for several even–even nuclei, again in both two-neutrino and neutrinoless modes. The lepton phase space for these transitions is suppressed compared to the decays of the type (1), while the nuclear structure issues are analogous. Hence, again we will concentrate on the decays described in equation (1). However, the 0νECEC mode in few cases can have a resonance character (Q → 0). The possible enhancement in those special cases is described in the appendix

The ββ decay, in either mode, can proceed only if the initial nucleus is stable against the standard β decay (both β⁻ and β⁺ or EC). That happens exclusively in even–even nuclei, where moreover the ground state is always I^π = 0⁺. The ββ decay rate is a steep function of the energy carried by the outgoing leptons (i.e. of the decay Q-value). Hence, transitions with larger Q-value are easier to observe. For this reason in table 1 I list all candidate nuclei with Q values larger than 2 MeV that are particularly well suited for the study of the ββ decay.

There has been a significant progress recently in the accuracy of the atomic mass determination using various trap arrangements. In many cases the Q-values are determined with accuracy better than 1 keV, making the search for the all important 0νββ decay mode easier; in table 1 these more recent Q-value determinations are shown, together with the corresponding references.

In both modes of the ββ decay the rate can be expressed as a product of independent factors that depend on the atomic physics (the so-called phase-space factors G^0ν and G^2ν) that include also the Q-value dependence as well as the fundamental physics constants, nuclear structure (the nuclear matrix elements M^0ν and M^2ν), and for the 0νββ mode the possible particle physics parameters (the effective neutrino mass 〈m_ββ〉 in the simplest case). Thus

$$\begin{equation} \frac{1}{T_{1/2}^{0\nu }} = G^{0\nu } |M^{0\nu }|^2 |\langle m_{\beta \beta } \rangle |^2 ; \hspace{28.45274pt} \frac{1}{T_{1/2}^{2\nu }} = G^{2\nu } |M^{2\nu }|^2 . \end{equation} \tag{ 2 }$$

The values of G^0ν and G^2ν are also listed in table 1. The entries there are taken from [6], and are not corrected for the small changes in Q and g_A since that time. Also, since by convention the nuclear matrix elements M^0ν are dimensionless, the nuclear radius appears in them as a multiplicative factor. To compensate for it, the phase-space factor G^0ν is proportional to R⁻², where R = r₀ × A^1/3 is the nuclear radius. In table 1 the value r₀ = 1.2 fm was used. (Note that, obviously, the values of the phase-space factors depend on the convention used for r₀ and g_A. One has to keep that issue in mind when using the equation (2) to relate the half-lives and nuclear matrix elements (see e.g. [15, 16]).)

Double beta transitions are possible and potentially observable because nuclei with even Z and N are more bound than the odd–odd nuclei with the same A = N + Z. A typical example is shown in figure 1. With one exception, all nuclei in table 1 have an analogous mass pattern. The one exception is ⁴⁸Ca where the intermediate nucleus ⁴⁸Sc can be in principle reached by the β⁻ decay of ⁴⁸Ca with Q = 278 keV. However, the ground state of ⁴⁸Sc is 6⁺ and the first excited state at 131 keV is 5⁺. β decays with a large nuclear spin change are heavily suppressed; in this particular case the β⁻ decay of ⁴⁸Ca has not been observed as yet, while the 2νββ decay has been observed.

Zoom In Zoom Out Reset image size

The two-neutrino mode (2νββ) is just an ordinary beta decay of two bound neutrons occurring simultaneously since the sequential decays are forbidden by the energy conservation law. For this mode, clearly, the lepton number is conserved and this decay is allowed in the standard model of electroweak interaction. It has been repeatedly observed in a number of cases and proceeds with a typical half-life of ∼10^19–20years for the nuclei with Q-values above 2 MeV. In contrast, the neutrinoless mode (0νββ) obviously violates the law of lepton number conservation and is forbidden in the standard model. Hence, its observation would, as already stated, be a signal of the 'physics beyond the standard model'.

The two modes of the ββ decay have some common and some distinct features. The common features are:

• The leptons carry essentially all available energy. The nuclear recoil is negligible, Q/Am_p ≪ 1.
• The transition involves the 0⁺ ground state of the initial nucleus and (in almost all cases) the 0⁺ ground state of the final nucleus. In few cases the transition to an excited 0⁺ or 2⁺ state in the final nucleus is energetically possible, but suppressed by the smaller phase space available. (But the 2νββ decay to the excited 0⁺ state has been observed in few cases.)
• Both processes are of second order of weak interactions, ∼G⁴_F, hence inherently slow. The phase space consideration alone (for the 2νββ mode ∼Q¹¹ and for the 0νββ mode ∼Q⁵) give preference to the 0νββ which is, however, forbidden by the lepton number conservation and therefore much slower (at least by a factor 10⁵ for the nuclei listed in table 1) than the 2νββ decay.

The distinct features are:

• In the 2νββ mode the two neutrons undergoing the transition are uncorrelated (but decay simultaneously) while in the 0νββ the two neutrons are correlated. In the 2νββ mode the corresponding momentum transfer q is restricted by the decay Q-value; hence qR ≪ 1. On the other hand, in the 0νββ mode the momentum transfer is of the order of the nucleon Fermi momentum q ∼ q_Fermi; hence qR ⩾ 1 in that case.
• In the 2νββ mode the sum electron kinetic energy T₁ + T₂ spectrum is continuous and peaked below Q/2. This is due to the electron masses and the Coulomb attraction. As T₁ + T₂ → Q the spectrum approaches zero approximately like (ΔE/Q)⁶.
• On the other hand in the 0νββ mode the sum of the electron kinetic energies is fixed, T₁ + T₂ = Q, smeared only by the detector resolution. This allows one to separate the two modes experimentally by measuring the sum energy of the emitted electrons with a good energy resolution, even if the decay rate for the 0νββ mode is much smaller than for the 2νββ mode.

Another hypothetical mode of double-beta decay is often considered in the literature, the decay accompanied by Majoron emission. Majorons are supposed to be very light or massless particles χ⁰ that couple to neutrinos. A variety of approaches involving massless (or almost massless) scalar particles and their impact on 0νββ decay have been considered. The discussion of this topic goes beyond the scope of the present review; a rather complete list of references can be found e.g. in [5]. In all of these hypotheses the sum electron spectra are continuous with the characteristic shape (neglecting Coulomb effects)

$$\begin{equation} \frac{{\rm d}\Gamma ^{0\nu M}}{{\rm d}E_1 {\rm d}E_2} \sim (Q - E_{1}\hbox{--}E_2)^n p_1 p_2 E_1 E_2 , \end{equation} \tag{ 3 }$$

where the 'index' n = 1–7 depends on the model in question. The decay rate is expressed as

$$\begin{equation} \Gamma ^{0\nu M} = | \langle g_{\chi } \rangle |^2 |M_{\chi }|^2 G_{\chi } (Q, Z) , \end{equation} \tag{ 4 }$$

where 〈g_χ〉 is the model dependent averaged coupling constant. For n = 1 that constant is experimentally constrained to be less than about 10⁻⁵.

The relation between the 0νββ-decay rate and the effective Majorana mass 〈m_ββ〉 is to some extent problematic. The rather conservative assumption leading to equation (2) is that the only possible way the 0νββ decay can occur is through the exchange of a virtual light, but massive, Majorana neutrino between the two nucleons undergoing the transition, and that these neutrinos interact by the standard left-handed weak currents. But that is not the only theoretically possible mechanism. Lepton number violating (LNV) interactions involving so far unobserved much heavier (∼ TeV) particles might lead to a comparable 0νββ decay rate.

In general 0νββ decay can be generated by (i) light massive Majorana neutrino exchange or (ii) heavy particle exchange (see, e.g. [17, 18]), resulting from LNV dynamics at some scale Λ above the electroweak one. The relative size of heavy (A_H) versus light particle (A_L) exchange contributions to the decay amplitude can be crudely estimated as follows [19]:

$$\begin{equation} A_L \sim G_F^2 \frac{\langle m_{\beta \beta } \rangle }{\langle k^2 \rangle } ,\quad A_H \sim G_F^2 \frac{M_W^4}{\Lambda ^5} , \quad \frac{A_H}{A_L} \sim \frac{M_W^4 \langle k^2 \rangle }{\Lambda ^5 \langle m_{\beta \beta } \rangle } \ , \end{equation} \tag{ 5 }$$

where 〈m_ββ〉 is the effective neutrino Majorana mass, 〈k²〉 ∼ (100 MeV)² is the typical light neutrino virtuality, and Λ is the heavy scale relevant to the LNV dynamics. Therefore, A_H/A_L ∼ O(1) for 〈m_ββ〉 ∼ 0.1–0.5 eV and Λ ∼ 1 TeV, and thus the LNV dynamics at the TeV scale leads to similar 0νββ-decay rate as the exchange of light Majorana neutrinos with the effective mass 〈m_ββ〉 ∼ 0.1–0.5 eV.

Obviously, the 0νββ lifetime measurement by itself does not provide the means for determining the underlying mechanism. The spin-flip and non-flip exchange can be, at least in principle, distinguished by the measurement of the single-electron spectra or polarization. However, in most cases the mechanism of light Majorana neutrino exchange, and of heavy particle exchange cannot be separated by the observation of the emitted electrons. We will not discuss here the possible ways of determining which mechanism is responsible for the transition, once the 0νββ decay has been observed. However, obviously, the corresponding nuclear matrix elements M^0ν depend on that.

Let us comment now on the main differences between the mechanism involving light or heavy particle exchange. We will show later that the evaluation of the nuclear matrix element M^0ν can be performed in the closure approximation, i.e. without explicit treatment of the virtual states in the intermediate odd–odd nucleus. Thus

$$\begin{equation} M^{0\nu } = \langle f || O^K || i \rangle , \end{equation} \tag{ 6 }$$

where the operator O^K creates two protons and annihilate two neutrons. In addition, the operator O^K depends on the distance between these nucleons, and on their other quantum numbers. The main difference between the mechanism involving the light massive Majorana neutrino exchange and the heavy (∼TeV) particle exchange is the range of the operator O^K. The light neutrino exchange represents two point-like vertices separated by the distance r ∼ 1/q. The decay rate is then proportional to the square of the effective Majorana neutrino mass as in equation (2). On the other hand the heavy particle exchange represents a single point-like vertex (six fermions, four hadrons and two leptons), i.e. dimension 9 operator. Proper treatment of the short range nucleon–nucleon repulsion is obviously crucial in that case. The relation between the neutrino mass and the decay rate is not simple in that case.

Independently of its mechanism the existence of 0νββ decay would mean that on the elementary particle level a six fermion LNV amplitude transforming two d quarks into two u quarks plus two electrons is nonvanishing. As was first pointed out by Schechter and Valle [20], already 30 years ago, this fact alone would guarantee that neutrinos are massive Majorana fermions. This qualitative statement (or theorem), however, as we pointed out above, does not in general allow one to deduce the magnitude of the neutrino mass once the rate of the 0νββ decay have been determined. It is important to stress, however, that quite generally an observation of any total LNV process, not only of the 0νββ decay, would necessarily imply that neutrinos are massive Majorana fermions.

Study of the 2νββ decay is an important nuclear physics problem by itself. Moreover, evaluation of the M^2ν matrix elements is an important test for the nuclear theory models that aim at the determination of the analogous but different quantities for the more fundamental 0ν neutrinoless mode. So, we begin our discussion with that experimentally more accessible ββ decay mode.

The rate of the 2νββ decay was first estimated by Maria Goeppert-Meyer already in 1937 in her thesis work suggested by E Wigner, basically correctly. Yet, first experimental observation in a laboratory experiment was achieved only in 1987, 50 years later [21]. (Earlier observations of the ββ decay [22–24] were based on the geochemical method that cannot separate the 2νββ and 0νββ decay modes. Later laboratory measurements have shown, as expected, that the geochemical method determines dominantly the 2νββ decay mode.)

Note that such delay is not really exceptional in neutrino physics. It took more than 20 years since the original suggestion of Pauli to show that neutrinos are real particles in the pioneering experiment by Raines and Cowan. And it took another almost 50 years since that time to show that neutrinos are massive fermions. Why it took so long in the case of the ββ decay? As pointed out above, the typical half-life of the 2νββ decay is ∼10^19–20 years. Yet, its 'signature' is very similar to natural radioactivity, present to some extent everywhere, and governed by the half-life of ∼10¹⁰ years, or much less for most of the man-made or cosmogenic radioactivities. So, background suppression is the main problem to overcome when one wants to study either of the ββ decay modes.

During the last two decades the 2νββ decay has been observed in 'live' laboratory experiments in many nuclei, often by different groups and using different methods. That shows not only the ingenuity of the experimentalists who were able to overcome the background nemesis, but makes it possible at the same time to extract the corresponding 2ν nuclear matrix element from the measured decay rate. The resulting nuclear matrix elements M^2ν, which have the dimension energy⁻¹, are plotted in figure 2. Note the pronounced shell dependence; the matrix element for ¹⁰⁰Mo is almost ten times larger than the ones for ¹³⁰Te or ¹³⁶Xe.

Zoom In Zoom Out Reset image size

The derivation of the 2νββ-decay rate formula is analogous to the treatment of ordinary beta decay. It begins with the Fermi golden rule for second order weak decay

$$\begin{equation} \frac{1}{\tau } = 2\pi \delta (E_0 - \Sigma _f E_f) \left[ \Sigma _{m,\mu } \frac{\langle f| H_{\mu } | m \rangle \langle m| H^{\mu } | i \rangle }{E_i - E_m - p_{\nu } - E_e} \right]^2, \end{equation} \tag{ 7 }$$

where the sum over m includes all relevant virtual states in the intermediate odd–odd nucleus and μ labels the different Dirac structures of the weak interaction Hamiltonian.

When calculating the 0⁺ → 0⁺ transitions, it is generally a very good approximation to replace the lepton energies in the denominator by the corresponding average value, i.e. E_e + p_ν ∼ E₀/2, where the E₀ = M_i − M_f is the total decay energy including electron masses. The lepton momenta, for both electrons and neutrinos are all q < Q and thus qR ≪ 1, where R is the nuclear radius. Hence the so-called long wavelength approximation is valid and the rate formula separates into a product of the nuclear and lepton parts, where the lepton part G^2ν(E₀, Z) is just the phase-space integral discussed and tabulated earlier.

The nuclear structure information is contained in the nuclear matrix element; only the Gamow–Teller στ part contributes in the long wavelength approximation

$$\begin{equation} M^{2\nu } = \Sigma _{m,i,k} \frac{ \langle 0^+_f | \vec{\sigma }_i \tau ^+_i | m \rangle \langle m | \vec{\sigma }_k \tau ^+_k | 0^+_i \rangle }{E_m - (M_i + M_f)/2} . \end{equation} \tag{ 8 }$$

The individual terms in the equation (8) have a well defined meaning, in particular for the most relevant ground state to ground state transitions. The terms $\langle m | \vec{\sigma }_k \tau ^+_k | 0^+_i \rangle$ represent the amplitudes of the β⁻ strength of the initial nucleus and can be explored in the nucleon charge exchange reactions such as (p, n) and (³He, t). On the other hand the terms $\langle 0^+_f | \vec{\sigma }_i \tau ^+_i | m \rangle$ represent the β⁺ strength in the final nucleus and can be explored in the opposite nucleon charge exchange reactions such as (n, p) and (d, ²He). In this way one can (up to the sign) explore the contribution of several low lying states to the M^2ν matrix element. (See the article by H Ejiri and D Frekers in this issue.)

In equation (8) the energy denominators signify that the low-lying intermediate I^π = 1⁺ states contribute significantly more that the high-lying states. Nevertheless, it is useful to analyze a related expression obtained in the closure approximation:

$$\begin{equation} M^{2\nu } = \Sigma _{m,i,k} \frac{ \langle 0^+_f | \vec{\sigma }_i \tau ^+_i | m \rangle \langle m | \vec{\sigma }_k \tau ^+_k | 0^+_i \rangle }{E_m - (M_i + M_f)/2} \equiv \frac{ \langle 0^+_f | \Sigma _{i,k} \vec{\sigma }_i \tau ^+_i \cdot \vec{\sigma }_k \tau ^+_k | 0^+_i \rangle }{\Delta \bar{E}} , \end{equation} \tag{ 9 }$$

where $\Delta \bar{E}$ is the average energy denominator defined by the above equation. The numerator of equation (9) is the definition of the closure 2ν nuclear matrix element. By itself, the 2ν closure matrix element does not have a direct physics interpretation.

For the two-body operator (i.e. in closure) in the 2ν matrix element one can derive an approximate sum rule, analogous to the famous Ikeda sum rule of the GT operator [25]. Remember that for the beta strengths $S^{\beta \mp } = \Sigma \vec{\sigma }_i \tau _i^{\pm }$ the Ikeda sum rule is obeyed independently of the nuclear structure as long as nuclei are assumed to be made of nucleons only

$$\begin{equation} S^{\beta -} - S^{\beta +} = 3(N-Z) . \end{equation} \tag{ 10 }$$

In nuclei with neutron excess, N > Z the β⁺ strength is much smaller than the β⁻ strength, which is concentrated in the giant GT resonance.

In analogy, the strength corresponding to the double beta operator $\Sigma _{i,k} \vec{\sigma }_i \tau _i^+ \vec{\sigma }_k \tau _k^+$ will be concentrated mostly in the 'double GT resonance' and its strength will be approximately equal to the sum rule 6(N − Z)(N − Z + 1) [25]. Assuming further that the average energy denominator $\Delta \bar{E}$ in equation (9) is O(1) MeV, we see that 2νββ closure matrix element, connecting the ground states of the initial and final nuclei (see figure 2) is heavily suppressed, representing only a very small fraction of the corresponding sum rule. From that it follows that an accurate evaluation of the matrix elements M^2ν and M^2ν_cl is going to be quite difficult, since these quantities are so small in their natural units.

To shed more light on the problem, consider the dependence of the 2νββ matrix elements on the distance between the two neutrons that are transformed into two protons. These two neutrons are not correlated, but nevertheless they are both bound in the corresponding nuclei, and decay together. To characterize the radial dependence let us introduce the function

$$\begin{eqnarray} &&C^{2\nu }_{{\rm cl}}(r) = \langle f | \Sigma _{lk} \vec{\sigma }_l \cdot \vec{\sigma }_k \delta (r - r_{lk}) \tau _l^+ \tau _k^+ | i \rangle , \nonumber \\ &&M^{2\nu }_{{\rm cl}} = \int _0^{\infty } C^{2\nu }_{{\rm cl}}(r)\, {\rm d}r . \end{eqnarray} \tag{ 11 }$$

This definition is in analogy to the related function C^0ν(r) first introduced in [26], and discussed in detail later. Note that while the matrix elements M^2ν and M^2ν_cl get contributions only from the 1⁺ intermediate states, the function C^2ν_cl gets contributions from all intermediate multipoles. This is the consequence of the δ function in the definition of C^2ν_cl(r). When expanded, all multipoles contribute. Naturally, when integrated over r only the contributions from the 1⁺ are nonvanishing. Examples of the function C^2ν_cl(r) are shown in figure 3. (Similar figure appears in [27].)

Zoom In Zoom Out Reset image size

One can see that in all the cases shown in figure 3 the function C^2ν_cl(r) consists of a positive peak at r ∼ 1 fm and a negative tail starting at r ∼ 2 fm. It turns out that the areas under the positive peak and the negative tail are nearly equal, resulting in considerable uncertainty in M^2ν_cl.

Theoretical evaluation of the matrix elements M^2ν and M^2ν_cl, respectively application of the corresponding rate equation (2), involves another problem, namely whether one should, or should not, apply the experience with the ordinary β decay and use the concept of quenching of the GT strength [28]. We shall discuss that issue in more detail later.

The 0ν decay rate associated with the nonvanishing value of m_ν is of the general form

$$\begin{equation} \omega _{0\nu } = 2 \pi \Sigma _{\rm spin} | R_{0\nu } |^2 \delta (E_{e1} + E_{e2} + E_f - M_i )\, {\rm d}^3 p_{e1} \,{\rm d}^3 p_{e2} , \end{equation} \tag{ 12 }$$

where E_f is the energy of the final nucleus and R_0ν is the transition amplitude including both the lepton and nuclear parts.

After substitution for the neutrino propagator and integration over the virtual neutrino momentum, the lepton amplitude acquires the form

$$\begin{equation} \fl {-}{\rm i} \delta _{jk} \int \frac{{\rm d}^4 q}{ (2 \pi )^4 } \frac{{\rm e}^{-{\rm i}q(x-y) }}{ q^2 - m_j^2 } \bar{e} (x) \gamma _{\rho } \frac{1}{2} (1 - \gamma _5 ) ( q^{\mu } \gamma _{\mu } + m_j ) \frac{1}{2} (1 - \gamma _5 ) \gamma _{\sigma } e^C (y) . \end{equation} \tag{ 13 }$$

From the commutation properties of the gamma matrices it then follows that the decay amplitude for purely left-handed lepton currents is proportional to the neutrino Majorana mass m_j. Integration over the virtual neutrino energy leads to the replacement of the propagator (q² − m²_j)⁻¹ by the residue π/ω_j with $\omega _ j = \big( \vec{q}^{ 2} + m_j^2 \big)^{1/2}$.

Finally, the integration over the space part ${\rm d} \vec{q}$ leads to an expression for the 'neutrino potential' that appears in the corresponding nuclear transition operator,

$$\begin{equation} \fl H(r,E_m ) = \frac{R}{ 2 \pi ^2} \int \frac{ {\rm d} \vec{q}}{ \omega } \frac{1}{ \omega + A_m } {\rm e}^{{\rm i} \vec{q} \cdot \vec{r}} = \frac{2R}{ \pi r} \int _0^ {\infty } {\rm d}q \frac{ q \sin (qr) }{ \omega ( \omega + A_m ) } = \frac{2R}{\pi } \int _0^ {\infty } {\rm d}q \frac{j_0(qr)q}{q + A_m} , \end{equation} \tag{ 14 }$$

where the nuclear radius R = 1.2A^1/3 fm was added as an auxiliary factor so that H becomes dimensionless. A corresponding 1/R² compensates for this auxiliary quantity in the phase space formula. The weak dependence on the excitation energy of the virtual intermediate odd–odd nucleus appears in A_m = E_m − E_i + E_e ≡ E_m − (M_i − M_f)/2.

The momentum of the virtual neutrino is determined by the uncertainty relation q ∼ 1/r, where r ⩽ R is a typical spacing between two nucleons. We will show later that in fact the relevant values of r are only r ⩽ 2–3 fm, so that the momentum transfer q ∼ 100–200 MeV. For the light neutrinos the neutrino mass m_j can then be safely neglected in the potential H(r). (Obviously, for heavy neutrinos, with masses M_j ≫ 1 GeV a different procedure is necessary.) Also, given the large value of q the dependence on the difference of nuclear energies E_m − E_i is expected to be rather weak and the summation of the intermediate states can be performed in closure for convenience. This approximation $H(r,E_m) \simeq H(r, \bar{E})$ is, in fact, typically used in the evaluation of the M^0ν, where $\bar{E} \sim 5\hbox{--}10$ MeV is the characteristic nuclear excitation energy.

It is worthwhile to test the validity of this approximation. Such test can be conveniently performed within the quasiparticle random phase approximation (QRPA), where the sum over the intermediate states can be easily explicitly carried out. In this context one can ask two questions: How good is the closure approximation? And what is the value of the corresponding average energy? In figure 4 we illustrate the answers to these questions (see also [27]). The QRPA matrix elements evaluated by explicitly summing over the virtual intermediate states quoted in the caption can be compared with the curves obtained by replacing all intermediate energies with a constant $\bar{E}$, which is varied there between 0 and 12 MeV. One can see, first of all, that the M^0ν changes modestly, by less than 10% when $\bar{E}$ is varied as expected given the relative sizes of q and $\bar{E}$ and, at the same time, that the exact results are quite close, but somewhat larger, than the closure ones. Thus, employing the closure approximation is appropriate for the evaluation of M^0ν even though it apparently slightly underestimates the M^0ν values.

Zoom In Zoom Out Reset image size

The neutrino potential in the equation (14) was defined assuming that the nucleons are point particles. That is not true, however, and thus it is necessary to include a corresponding correction in the definition of $H(r,\bar{E})$. It is customary to approximate this correction in the form of the dipole type form factor

$$\begin{equation} f_{\rm FNS} = \frac{1}{\big(1 + \frac{q^2}{M_A^2}\big)^2} , \end{equation} \tag{ 15 }$$

with M_A = 1.09 GeV. Varying M_A in the interval 1.0–1.2 GeV makes little difference.

In addition, while the weak current of quarks has the simple V–A structure, the weak nucleon current contains additional terms, since nucleons are complicated composite objects,

$$\begin{equation} \fl J^{{\mu }^+} = \Psi \tau ^+ \left[ g_V(q^2)\gamma ^{\mu } - {\rm i}g_M(q^2)\frac{\sigma ^{\mu \nu }}{2 m_p} q_{\nu } -g_A(q^2)\gamma ^{\mu }\gamma _5 + g_P(q^2)q^{\mu }\gamma _5 \right] \Psi , \end{equation} \tag{ 16 }$$

representing the vector, weak magnetism, axial vector and induced pseudoscalar. In the zero momentum transfer limit the values of the corresponding form factors are well known. The induced pseudoscalar form factor is usually evaluated using the partially conserved axial-vector current hypothesis

$$\begin{equation} g_P(q^2) = \frac{g_A(q^2) \times 2m_p}{q^2 + m_{\pi }^2} , \end{equation} \tag{ 17 }$$

and the weak magnetism form factor is simply proportional to g_V(q²), g_M(q²) = (μ_p − μ_n)g_V(q²).

Taking these 'higher order terms' into account and using the nonrelativistic limit for the nucleons, one arrives at the corresponding correction term g_HOT(q²) in the main, GT part of the neutrino potential (for details see [29])

$$\begin{equation} \fl g_{\rm HOT}^{{\rm GT}}(q^2) = 1 - \frac{2}{3}\frac{\vec{q}^2}{\vec{q}^2 + m_{\pi }^2} + \frac{1}{3} \left( \frac{\vec{q}^2}{\vec{q}^2 + m_{\pi }^2} \right)^2 + \frac{2}{3} \left(\frac{g_V(q^2)}{g_A(q)} \right)^2 \frac{(\mu _p - \mu _n)^2 \vec{q}^2}{4 m_p^2} . \end{equation} \tag{ 18 }$$

Finally, the neutrino potential $H_{{\rm GT}}(r,\bar{E})$ governing the Gamow–Teller part of the matrix element M^0ν with these correction factors included is of the form

$$\begin{equation} H_{{\rm GT}}(r,\overline{E}_{0\nu }) = \frac{2R}{\pi } \int _0^{\infty } j_0(q r)\frac{q}{q + \overline{E}_{0\nu }} f^2_{\rm FNS}(q^2)g_{\rm HOT}^{{\rm GT}}(q^2)\,{\rm d}q . \end{equation} \tag{ 19 }$$

When the finite nucleon size, higher order terms are neglected, and $\bar{E}_{0\nu } = 0$ is assumed, the potential has Coulomb-like shape R/r. The full potential, equation (19), however, is finite at r = 0, $H(r \rightarrow 0, \bar{E}_{0\nu }=0) = 5M_A R/16$. Including the higher order currents and finite $\bar{E}$ in equation (19) increases the value of H(r = 0) by ∼30%. The shape of the neutrino potential is shown in figure 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the full 0νββ operator there are several potentials, each with its own spin structure

$$\begin{equation} \fl O^{0\nu } = \Sigma _{i,k} \tau _i^+ \tau _k^+ \left(-H_F(r,\overline{E}_{0\nu })\frac{g_V^2}{g_A^2} + H_{{\rm GT}}(r,\overline{E}_{0\nu }) \vec{\sigma }_i \cdot \vec{\sigma }_k + H_{T}(r,\overline{E}_{0\nu })S_{ik} \right), \end{equation} \tag{ 20 }$$

where in the momentum representation $S_{ik} = 3 \vec{\sigma }_i \cdot \vec{q} \vec{\sigma }_k \cdot \vec{q} - \vec{\sigma }_i \cdot \vec{\sigma }_k$ and g^F_HOT(q²) = 1, $g_{\rm HOT}^{T}(q^2) = \frac{2}{3} \frac{\vec{q}^2}{\vec{q}^2 + m_{\pi }^2} - \frac{1}{3} \big( \frac{\vec{q}^2}{\vec{q}^2 + m_{\pi }^2} \big)^2$. Also, in f_FNS for the Fermi part one needs to replace M_A by M_V ≃ 0.85 GeV, and in the tensor part g^T_HOT(q²) the spherical Bessel function j₀(qr) must be replaced by j₂(qr).

To finish this section, let us briefly mention the analysis of the 0νββ decay assuming the existence of the right-handed weak currents. Such phenomenological approach is often quoted in the literature, even though it is not at all clear that there exist a corresponding realistic particle physics model giving the 0νββ-decay rates competitive with the standard light left-handed Majorana neutrino exchange, like in equation (2). The assumed Hamiltonian is

$$\begin{equation} H_W = \frac{G_F \cos \theta _C}{\sqrt{2}} \Sigma _{i=1}^{2n} \big[ j_{Li\mu } \bar{J}_{Li}^{\mu \dag } + j_{Ri\mu } \bar{J}_{Ri}^{\mu \dag } \big] , \end{equation} \tag{ 21 }$$

where the lepton currents are

$$\begin{equation} j_{Li}^{\mu } = \bar{e}(1 - \gamma _5)N_{iL} ;\qquad j_{Ri}^{\mu } = \bar{e}(1 + \gamma _5)N_{iR} , \end{equation} \tag{ 22 }$$

with the hypothetical heavy neutrinos N_iL(R) with n mass eigenstates. The nuclear currents are

$$\begin{equation} \bar{J}_{Li}^{\mu \dag } = U_{ei} J_L^{\mu \dag } ; \qquad\bar{J}_{Ri}^{\mu \dag } = V_{ei} \big[\lambda J_R^{\mu \dag } + \eta J_L^{\mu \dag } \big] , \end{equation} \tag{ 23 }$$

where U_ei is the standard light neutrino mixing matrix and V_ei is its analogue for the heavy neutrinos. Parameters λ and η characterize the strength of the right–right and right–left interactions. The formulae for the corresponding 0νββ-decay rate have been evaluated in many papers, here I quote just two of them [30, 31]. In those references one can find the complete expressions for all necessary neutrino potentials and the corresponding spin and $\vec{q}$ dependencies.

As stressed above, when heavy particles of any kind mediate the 0νββ decay, we are dealing with a six fermion vertex, representing extremely short range operator. Vergados [32] was presumably the first author to describe how the issue of suppression, due to the short range nucleon–nucleon repulsion, can be overcome.

If one could treat nucleons as pointlike particles and assume that the heavy neutrino mass will be substantially larger than the corresponding momentum transfer q (say, m_heavy ≫ 1 GeV) the neutrino potential would contain exp(−mr) where m = min(m_heavy, M_W). The corresponding nuclear matrix element would be heavily suppressed due to the extremely short range of this potential.

However, nucleons are not pointlike particles. Instead, their finite size could be parametrized through the dipole type form factor f_FNS in equation (15). In that case the neutrino potential will be of the form

$$\begin{equation} \fl H_{{\rm heavy}\nu }(r, M_A) = \frac{4\pi R}{M_A^2} \int \frac{{\rm d} \vec{q}}{(2\pi )^3} \left( \frac{M_A^2}{M_A^2 + \vec{q}^2} \right)^4 = \frac{M_A R}{48} {\rm e}^{-M_A r} \left[ 1 + M_A r + \frac{1}{3} (M_a r)^2 \right] . \end{equation} \tag{ 24 }$$

Such potential will have the range 1/M_A and will be much less affected by the short range correlations (see e.g. [32, 31]). The disadvantage of the form factor modeling is that the error introduced by such approximation is very difficult to estimate.

In general the heavy particle exchange is characterized by some energy scale Λ_ββ (see equation (5)), that is much larger than any hadronic scale Λ_H ∼ 1 GeV that would enter the problem. In that case one could, instead of the form factor approach indicated above, use the prescriptions of effective field theory and classify the contributions in powers of small quantities q/Λ_H, q/Λ_ββ and Λ_H/Λ_ββ as in [33].

The hadronic vertices appearing in the corresponding Lagrangian will be of the type NNNNee, NNπee and ππee as illustrated in figure 7. They stem from quark–lepton operators having different transformation properties under parity and chiral SU(2). As such, they will contribute to different orders in the q/Λ_H expansion. The vertices involving pions are longer range. They have been analyzed in the form factor approach in [34], but the EFT allows more systematic approach because of the separation of scales q < Λ_H ≪ Λ_ββ. It was noted already in [35] that the nuclear matrix elements associated with the long range pionic effects within the RPV SUSY scenarios can be dominant. But that is, in fact, a more general result. The pionic effects can be substantially larger than those obtained using the conventional form factor model for the short-range NNNNee process.

Zoom In Zoom Out Reset image size

The LNV vertices in figure 7 represent nonstandard model operators, and can be characterized by some parameters K_ππ, K_NNπ and K_NNNN. Counting powers of q/Λ_H, for Λ_H = 4πf_π, f_π ≃ 92 MeV, we come to the conclusion that

$$\begin{equation} \fl {\rm Figure 7(a)}\sim K_{\pi \pi } p^{-2},\qquad {\rm figure 7(b)} \sim K_{NN\pi }p^{-1},\qquad {\rm figure 7(c)} \sim K_{{NNNN}}p^0. \end{equation} \tag{ 25 }$$

Thus, the long range 0νββ-decay operators ππee and to a lesser degree NNπee are enhanced, relative to the short range operator NNNNee in figure 7(c).

Clearly, the most important operator is the one corresponding to the ππee vertex. It was shown in [33] that the corresponding Lagrangian can be expressed as

$$\begin{equation} {\mathcal {L}}_{(0)}^{\pi \pi ee} = \frac{G_{{\rm F}}^2\Lambda _{{\rm H}}^2 f_\pi ^2}{\Lambda _{\beta \beta }}\lbrace \pi ^-\pi ^-\bar{e}(\beta _1 + \beta _2\gamma ^5)e^c + \pi ^+\pi ^+\bar{e}^c(\beta _1 - \beta _2 \gamma ^5)e \rbrace , \end{equation} \tag{ 26 }$$

where β_{1, 2} are dimensionless parameters that need be evaluated in any concrete particle physics model.

Transforming the above Lagrangian into the hadron–lepton system, taking the nonrelativistic limit and Fourier transforming to coordinate space yield an expression for the transition matrix element

$$\begin{eqnarray} &&{\mathcal {M}}_0 = \langle\Psi _{A,Z+2}| \sum \limits _{ij}\frac{R}{r_{ij}} [ F_1\!(x_{ij}) \vec{\sigma }_i \cdot \vec{\sigma }_j + F_2\!(x_{ij}) S_{ij} ] \tau ^+_i\tau ^+_j | \Psi _{A,Z} \rangle, \\ &&\nonumber S_{ij} = 3\vec{\sigma }_i \cdot \hat{\vec{r}}_{ij} \vec{\sigma }_j \cdot \hat{\vec{r}}_{ij} - \vec{\sigma }_i \cdot \vec{\sigma }_j . \end{eqnarray} \tag{ 27 }$$

Here r_ij is the distance between the ith and jth neutrons, $\hat{\vec{r}} = \vec{r}/r$ and x_ij = r_ijm_π. The form-factors F₁ and F₂ were first introduced in [34]

$$\begin{equation} {\textit {F}}_1(x)=(x-2)\,{\rm e}^{-x},\qquad {\textit {F}}_2(x)=(x+1)\,{\rm e}^{-x} . \end{equation} \tag{ 28 }$$

The 0νββ half-life is, as expected, inversely proportional to the square of the scale Λ_ββ and depends on the empirical parameters β_{1, 2},

$$\begin{eqnarray} &&\fl \frac{1}{T_{1/2}} = \frac{\hbar c^2}{144\pi ^5\ln 2}\frac{g_{{\rm A}}^4}{R^2} \frac{\Lambda _{H}^4 G_F^4}{\Lambda _{\beta \beta }^{2}} \int _{m_e}^{E_{\beta \beta } - m_e}{\rm d}E_1 F\!(Z+2,E_1)F\!(Z+2,E_2) \nonumber \\ &&\times \frac{1}{2} \big[\big(\beta _1^2+\beta ^2_2\big)p_1E_1p_2E_2 - \big(\beta _1^2-\beta ^2_2\big)p_1p_2m_e^2\big] |{\mathcal {M}}_0|^2. \end{eqnarray} \tag{ 29 }$$

As mentioned earlier, concrete application of the pion exchange mechanism in [35] suggests the dominance of the ππee over the short range nucleon only vertex by a factor of 10–30 at the level of the nuclear matrix element. This is, to some extent, so far largely unexplored issue. The application of the form factor method, obviously insufficient, is more or less a norm even today (see e.g. the recent paper [36]).

When evaluating the 0νββ nuclear matrix element in the closure approximation it is necessary, ultimately, to consider the matrix element of a two-body transition operator which connects a state with two neutrons in a specific mean field single particle states with two protons again in some specific final states. These states are angular momentum and possibly isospin coupled single particle states. In essentially all cases they do not explicitly include the effect of the short range nucleon–nucleon repulsion.

Until recently, it was customary to simulate the effect of such short range correlations by including in the operator (or equivalently in the two body wave-function) a phenomenological Jastrow-like function

$$\begin{equation} f_{\rm Jastrow}(r) = 1 - c {\rm e}^{-ar^2}(1 - br^2) . \end{equation} \tag{ 30 }$$

In other words,

$$\begin{equation} M^{0\nu } = \Sigma _{k<m} M_{km} M_{km} \rightarrow f(r_{km})M_{km} f(r_{km}) . \end{equation} \tag{ 31 }$$

In early works the parameters a, b, c, determined for a different purpose in [37], were used (c = 1.0, a = 1.1 fm⁻² and b = 0.68 fm⁻²). This is the so-called Miller–Spencer parametrization. Since in the evaluation of the nuclear matrix element such function appears twice, we plot in figure 8 its square. Application of this prescription resulted typically in the reduction of the corresponding 0νββ nuclear matrix element by ∼30% when compared to the result when f_Jastrow(r) = 1 is used.

Zoom In Zoom Out Reset image size

Recent works have questioned the adequacy of this prescription. One of the approaches proposed relatively recently in [38] is based on the unitary correlation operator method (UCOM), see [39]. This approach describes not only short-range, but also central and tensor correlations explicitly by means of a unitary transformation. Applied to a realistic NN interaction, the method produces a correlated interaction, which can be used as a universal effective interaction. In the case of the 0νββ-decay calculation the correlated two-nucleon wave function was taken as

$$\begin{equation} | {{\Psi }}_{\mathcal J} \rangle = C_r | {{\Psi }}_{\mathcal J} \rangle . \end{equation} \tag{ 32 }$$

Here, C_r is the unitary correlation operator describing the short-range correlations. The explicit form of C_r is given in [39]. The UCOM-corrected nuclear matrix elements are significantly less reduced when compared with those calculated with Jastrow SRC [40, 38]. In practice, the UCOM prescription requires that the operator r_ab be replaced by a shifted R₊(r_ab). As pointed out in [41], care must be taken to apply the prescription consistently, not just simply by replacing the neutrino potential H(r) → H(R₊(r)).

In [42, 43] the effect of short range correlations was computed within well defined Brueckner-based approximation scheme. In particular, in [43] the coupled cluster method was used, based on the realistic NN interactions, CD-Bonn and Argonne V18. In a good approximation the effect of correlations can be again approximated by the Jastrow-like functions (see equation (30)) with parameters fitted in [43], a = 1.59 fm⁻², b = 1.45 fm⁻² and c = 0.92 for the case of the Argonne potential and a = 1.52 fm⁻², b = 1.88 fm⁻² and c = 0.46 for the case of the Argonne potential. Both of these functions are displayed in figure 8. It is obvious that the application of these functions leads to the much less reduction of the M^0ν, in agreement with the UCOM method and with the conclusion of [42] than the Miller–Spencer approach used previously. In fact, provided the nucleon finite size form factor are properly taken into account, the application of the UCOM procedure or the prescripition of [43] leads to the M^0ν values essentially unchanged in comparison to the omission of the short range correlation correction altogether. Thus, the M^0ν are about 30% larger than after the application of the traditional Miller–Spencer parametrization. There is a general consensus that the matter of treatment of the short range correlation is satisfactorily resolved at the present time.

In this section the basic issues involved in the numerical evaluation of the 0νββ nuclear matrix elements will be described. We concentrate here on the standard light Majorana neutrino exchange mechanism, so that the relation between the decay rate and the nuclear matrix elements is governed by equation (2). First, there are few reasonable assumptions common to all (or essentially all) methods:

• In the 0⁺ → 0⁺ transitions the outgoing leptons are in the s_1/2 (or more precisely Dirac's κ = −1) state. That assumption is used in the calculation of the phase space factor.
• The hadronic currents are treated in the nonrelativistic impulse approximation.
• In most, but not all, applications the closure approximation is used, i.e. the matrix element of a two-body operator is considered (see figure 4).

Since the decay involves the transformation of two bound neutrons into two bound protons, it is necessary, first of all, to choose the proper mean field in which the nucleons are bound. That field could be spherical, but it can be also deformed. In the case of deformed mean field the corresponding intrinsic states have no definite angular momentum; projection into states with definite value of J^π is sometimes performed in that case.

Once the mean field has been specified, the set of single particle states |ψ_i〉 is obtained. The next step is to decide which of these single particle states are fully occupied or totally empty in both the initial and final nuclei. Such states, obviously, cannot participate in the decay; the fully occupied states form the inert core. The remaining states, i.e. the single particle states that are partially occupied or where the occupancies in the initial and final nuclei are different, form the valence space. To obtain the resulting many-particle states one has to take into account the residual interaction among the nucleons constrained in the valence space.

While such division into occupied, valence, and empty states appears to be reasonable and well defined, in practice real nuclei have diffuse Fermi levels and it is not a priori clear how to properly decide to which category a given single particle state belongs. Moreover, the effect of the core and empty states should be included, in principle, in the renormalization of the nuclear Hamiltonian and in the definition of the effective operators. See [42] and the contribution of J Engel to this focus issue for explanation and references; there are not many explicit applications of this procedure to the problem at hand, so far. Ideally, moreover, the mean field should be determined selfconsistently, so that the same Hamiltonian, or nucleon–nucleon interaction, is used in its definition and in the treatment of the interaction of valence nucleons. Note that for ease of computation the single particle wave functions are typically taken to be the eigenfunctions of the harmonic oscillator potential, or superpositions of such functions.

Variety of methods has been used for evaluation of the 0νββ nuclear matrix elements. They differ in their choice of the valence space, interaction Hamiltonian and the ways the corresponding equations of motion are solved. Some of these methods are described in detail in separate contributions to this focus issue. Here I wish to characterize each method very briefly and show its typical output. Let me stress that an exact 'ab initio', i.e. without approximations, calculation of M^0ν for the candidate nuclei is impossible at the present time.

The nuclear shell model (NSM) is, in principle, the method that seems to be well suited for this task. In it, the valence space consists of just few single particle states near the Fermi level (usually one main shell). With interaction that is based on the realistic nucleon–nucleon force, but renormalized slightly to describe better masses, energies and transitions in real nuclei, all possible configurations of the valence nucleons are included in the calculation. The resulting states have not only the correct number of protons and neutrons, but also all relevant quantum numbers (angular momentum, isospin, etc). For most nuclei of interest (⁴⁸Ca is an exception) the valence space, however, does not include enough states to fulfil the Ikeda sum rule (see equation (10)), hence full description of the β strength functions S^{β ±} is not possible. However, NSM is well tested, since it is capable to describe quite well the spectroscopy of low lying states in both initial and final nuclei. In the following figures 9–11 the NSM results are denoted by the blue squares.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The 2νββ decay matrix elements M^2ν for several nuclei in table 1 are reasonably well described in the NSM, see [44] (¹⁰⁰Mo being a notable exception). However, to achieve this task, it was necessary to apply quenching factors that, for nuclei heavier than ⁴⁸Ca, are considerably smaller than in the lighter nuclei where the valence space contains the full oscillator shell. Note that no quenching is applied to the results shown in figures 9–11. I will describe the issue of quenching of the weak nucleon current operators in section 9.

The QRPA and its renormalized version (RQRPA) is another method often used in the evaluation of M^0ν. In it, the valence space is not restricted and contains at least two full oscillator shells, often more than that. On the other hand, only selected simple configurations of the valence nucleons are used. The basis states have broken symmetries in which particle numbers, isospin, and possibly angular momentum are not good quantum numbers but conserved only on average. After the equations of motion are solved, some of the symmetries are partially restored. The RQRPA partially restores the Pauli principle violation in the resulting states.

The procedure consist of several steps. In the first one the like particle pairing interaction is taken into account, using the BCS procedure. Then, the neutron–proton interaction is used in the equations of motion, resulting in states that contain two quasiparticle and two quasihole configurations and their iterations. Usually, the realistic G-matrix based interaction is used, but block renormalized with common renormalization factors. In particular, all particle–particle channel interaction matrix elements are typically scaled by the factor g_pp with nominal, unrenormalized value g_pp = 1. In many applications, including the one used in figures 9–11, the parameter g_pp is determined such that the experimental matrix element M^2ν is correctly reproduced. This procedure has been first suggested in [45], resulting in <20% change in g_pp. In this case, the method obviously is not predicting the magnitude of M^2ν, instead it uses its value in the fit. However, as shown in [45], the procedure removes the dependence of M^0ν on the size of the valence space, making its predicted value essentially constant. While most early QRPA and RQRPA calculations assume spherical nuclear shape, recent [46] extension allows deformed nuclear shape. In figures 9–11. the RQRPA results are denoted by red circles.

The IBM-2 method uses the microscopic interacting boson model to evaluate M^0ν [47]. In IBM-2 one begins with correlated S and D pairs of identical nucleons and includes the effect of deformation through the bosonic neutron–proton quadrupole interaction. The method describes well the low lying states, the electromagnetic transitions between them and the two-nucleon transition rates in spherical and strongly deformed nuclei. Even though the method was originally considered as an approximation of the NSM, the resulting M^0ν (purple upward pointing triangles in figures 9–11) are rather surprisingly close to the RQRPA results and noticeably larger than the NSM ones. At the present time it is not possible to evaluate the M^2ν matrix elements within the IBM-2 method.

The projected Hartree–Fock–Bogoljubov framework [48] (PHFB) uses angular momentum projected wave functions based on the Hamiltonian with pairing and quadrupole–quadrupole interaction. Again, the method is not capable in describing the intermediate odd–odd nuclei and thus neither the M^2ν matrix elements. The resulting M^0ν matrix elements for the heavier nuclei (including the deformed ¹⁵⁰Nd) are denoted with brown downward pointing triangles in figures 9–11. This method, like IBM-2, does not explicitly contain the isoscalar neutron–proton interaction, crucially important in QRPA. Yet its resulting M^0ν are in a reasonable agreement with those from RQRPA.

In [49] the generating coordinate method (GCM) was used in conjunction with the particle number and angular momentum projection for both the initial and final nuclei. Large single particle space was used (11 shells) with the well established Gogny D1S energy density functional. The initial and final many-body wave functions are represented as combinations of the particle number N, Z projected I = 0⁺ axially symmetric states with different intrinsic deformations. The lowest state is found by solving the Hill–Wheeler–Griffin eigenvalue equation. The resulting M^0ν matrix elements are denoted with green diamonds in figures 9–11. Again, like in PHFB and IBM-2 methods, the treatment of the odd–odd nuclei, and thus also of the M^2ν matrix elements is impossible in GCM.

Looking at the figure 9 one can see that, common to all five displayed methods the predicted M^0ν nuclear matrix elements vary relatively smoothly, with the mass number A, unlike the experimentally determined M^2ν matrix elements displayed in figure 2 that vary strongly with A. The RQRPA, IBM-2, PHFB and GCM method are in a crude agreement with each other, and predict slow decrease of M^0ν with increasing A. The M^0ν evaluated in NSM are essentially constant with A and noticeably smaller than those from the other methods, particularly in the lighter nuclei. In figures 10 and 11 I show the half-lives T_1/2 for the fixed effective neutrino Majorana mass 〈m_ββ〉 and 〈m_ββ〉 for fixed T_1/2. These figures should be helpful for comparing the experimental results and/or sensitivities for different nuclei.

The question that is, for obvious reasons, often asked is what is the error or uncertainty of individual nuclear matrix elements? It is notoriously difficult to estimate the error of a theoretical result. One possibility is just to use, without selection, all calculated values and take as the error their spread. That was used in the provocative [50]. Clearly that is not the correct approach, even though the quoted paper served as appropriate warning to the nuclear structure community.

Another possibility is to use as a measure of uncertainty the spread of selected careful calculations, like those in figure 9. That also has problems. Different methods use different approximations, and it is unlikely that all of them are equally important. And it is difficult to decide which of them is more realistic and which is less realistic.

Within a given nuclear structure method it is possible to assign a measure of error by considering the uncertainties in the input parameters. That was done for QRPA e.g. in [45, 40] where the size of the single particle basis, whether QRPA or RQRPA was used, and whether g_A was quenched or not is used for such estimate. Clearly that error does not include the possible systematic uncertainty of the basic method itself, but serves as a useful estimate of the irreducible uncertainty. When considering the ratios of the matrix elements for different nuclei, the question of correlation of the corresponding errors arises. As shown in [51] at least in the case of QRPA, the errors are highly correlated.

Double beta decay in both 2ν and 0ν modes can exist because even–even nuclei are more bound than the neighboring odd–odd nuclei. This extra binding is a consequence of pairing between like nucleons. In nonmagic systems neutrons and/or protons form 0⁺ pairs and the corresponding Fermi level becomes diffuse over the region with the characteristic size ∼ pairing gap Δ. This opens more possibilities for nn → pp transitions. The calculated matrix elements M^0ν increase when the gap Δ increases. The (unrealistic in real nuclei) situation of pure paired nuclear system (only seniority 0 states) would have very large M^0ν.

However, in real nuclei opposite tendencies exist. Real nuclei have admixtures of the 'broken pair' states, or in the shell model language, states with higher seniority. These states are present because other parts of the nucleon interaction exist, in particular the neutron–proton force. It is illustrative to characterize such states by the angular momentum $\mathcal {J}$ of the neutron pair that is in the ββ decay transformed into the proton pair with the same $\mathcal {J}$. In figure 12 the corresponding competition is illustrated, in both NSM and QRPA. While the positive pairing parts are very similar (since the same single particle spaces are used), the negative $\mathcal {J} \ne$ 0 parts are only qualitatively similar. Nevertheless the severe cancellation between these two tendencies, and between the corresponding components of the residual interaction, is present in both methods.

Zoom In Zoom Out Reset image size

The situation shown in figure 12 is somewhat unrealistic as far as the QRPA is concerned; the corresponding single particle space is just too small. When it is enlarged to contain two full oscillator shells, the competition becomes even more severe. Both the $\mathcal {J}$ = 0 and $\mathcal {J} \ne$ 0 grow, while their difference remains more or less unchanged. This situation, illustrated in figure 13, appears to be is universal. One can see that, due to this cancellation, accurate evaluation of M^0ν is difficult; small uncertainty in either of the two (positive and negative) components is strongly enhanced in the difference.

Zoom In Zoom Out Reset image size

The 0νββ decay operators depend on the internucleon distance r₁₂ due to the neutrino potential $H(r,\bar{E})$. Obviously, the range of r₁₂ is restricted from above by r₁₂ ⩽ 2R. From the form of H(r) ∼ R/r one could, naively, expect that the characteristic value of r₁₂ is the typical distance between nucleons in the nucleus, namely $\bar{r}_{12} \sim R$. However, that is not true, in reality only much smaller values of r₁₂ ⩽ 2–3 fm or equivalently larger values of the momentum transfer q are relevant. That was first demonstrated within QRPA in [40], very similar result was obtained also in NSM [52].

To see how that conclusion is obtained, define a function C(r)

$$\begin{equation} C^{0\nu }_{{\rm GT}}(r) = \langle f | \Sigma _{lk} \vec{\sigma }_l \cdot \vec{\sigma }_k \tau _l^+ \tau _k^+ \delta (r - r_{lk}) H(r_{lk},\bar{E}) | i \rangle . \end{equation} \tag{ 33 }$$

Obviously, this function is normalized by

$$\begin{equation} M^{0\nu }_{{\rm GT}} = \int _0^{\infty } C^{0\nu }_{{\rm GT}} (r)\, {\rm d}r , \end{equation} \tag{ 34 }$$

and is related to the function C^2ν_cl(r) (see equation (11)) by

$$\begin{equation} C^{0\nu }_{{\rm GT}}(r) = H(r,\bar{E}) \times C^{2\nu }_{{\rm cl}}(r) , \end{equation} \tag{ 35 }$$

which is valid for any shape of the neutrino potential.

Examples of the functions C^0ν(r) are shown in figure 14 for three representative nuclei. As the lower panel demonstrates, the cancellation between the 'pairing' ($\mathcal {J}$ = 0) and 'broken pairs' ($\mathcal {J} \ne$ 0) is essentially complete for r₁₂ ⩾ 2–3 fm (similar figure appears in [40]). That cancellation is so complete when the particle–particle renormalization constant g_pp is chosen in the usual way, i.e. so that the 2νββ lifetime is correctly reproduced. For other values of g_pp the cancellation between $\mathcal {J}$ = 0 and $\mathcal {J} \ne$ 0 is less perfect. Analogous conclusions can be obtained in an exactly solvable model [26] based on the algebra SO(5) × SO(5).

Zoom In Zoom Out Reset image size

The relation between C^0ν(r) and C^2ν_cl(r) in equation (35) is intriguing. Can it be used for the evaluation of M^0ν? There are two problems with that possibility. While M^2ν can be deduced from the measured lifetime of the 2νββ decay, the closure matrix element cannot be deduced that way. Moreover, even if it would be possible to somehow determine M^2ν_cl accurately, its magnitude is insufficient to determine the whole function C^2ν_cl(r).

It is well known that experimental Gamow–Teller β-decay transitions to individual final states are noticeably weaker than the theory predicts. That phenomenon is known as the axial current matrix elements quenching. The β-strength functions can be studied also with the charge exchange nuclear reactions and a similar effect is observed as well. Thus, in order to describe the matrix elements of the operator στ, the empirical rule (στ)²_eff ≃ 0.6(στ)²_model is usually used (see [53–55]). Since these operators accompany weak axial current, it is convenient to account for such quenching by using an effective coupling constant g^eff_A ∼ 1.0 instead of the true value g_A = 1.269.

The evidence for quenching is restricted so far to the Gamow–Teller operator στ and relatively low-lying final states. It is not known whether the other multipole operators associated with the weak axial current should be quenched as well. In fact, the analysis of the muon capture rates in [56] and of the unique second forbidden β decays in [57] suggests that quenching is not needed in those cases.

The quenching of the axial current matrix elements, or simply of g_A, is believed to be caused by nucleon correlations [58]. The alternative explanation, screening of the GT operator by the Δ-hole pairs [59] is, these days, not considered as the likely explanation. Very recently, in [60], the quenching has been associated with the two-body currents.

Since the 2νββ decay involves only the GT operators and relatively low-lying intermediate states, one could expect that the quenching is involved in that case. Indeed, as already mentioned, in the NSM the agreement with experimental decay rate is achieved only with quenching (στ)²_eff ≃ (0.2 − 0.36)(στ)²_model, with similar quenching required to reproduce the measured β decay rates and the β strength functions.

In QRPA and RQRPA the usual way of the g_pp renormalization using the experimental M^2ν means that the effect of quenching cannot be deduced from the comparison of the calculated and empirical M^2ν. In the other discussed methods (IBM-2, PHFB, GCM) it is not possible at the present time to calculate M^2ν, so the question of quenching for the 2νββ is meaningless in that case.

However, whether quenching should be applied to the 0νββ matrix elements M^0ν is an important issue and a source of noticeable uncertainty. Within QRPA and RQRPA the effect of quenching on the 0νββ decay mode has been considered in [27, 43, 61]. The calculations there uses a modified definition

$$\begin{equation} M^{\prime 0\nu } = \left(\frac{g^{\rm eff}_A}{1.269}\right)^2 M^{0\nu }\big(g^{\rm eff}_A\big), \end{equation} \tag{ 36 }$$

and thus the standard phase space factor can be used independently of whether quenching is included or not. It is concluded there that the matrix elements M'^0ν are reduced by ∼20–30% when g^eff_A = 1.0 is assumed compared with g_A = 1.269. Thus, in that case the predicted 0ν decay rates are affected by the possible quenching less than naive expectation based on the ratio [g^eff_A/g_A]⁴ might suggest.

Somewhat similar, but from the point of view of physics involved quite different reasons was reached in [60] when the effect of two-body current is included in the effective field theory.

The problem of quenching of the M^0ν matrix elements remains one of the main sources of uncertainty in their true value. The relative change with A and Z, as depicted in figure 9 will be, however, affected only little by the inclusion of the quenching phenomenon.

In the last few years the problem of evaluation of the nuclear matrix elements M^0ν of 0νββ decay has entered a new era. As described here, there are now five or more seemingly quite different methods of solving it; it became one of the very active subfields of nuclear structure theory. The results are not yet quite convergent, but at least they suggest that the calculated values are fairly insensitive (within a factor of about 2) to the broad range of approximations made. Moreover, all calculations agree that the values of M^0ν do not change abruptly from one candidate nucleus to another one. Thus, if the 0νββ decay were observed in one nucleus, one can with some confidence predict its lifetime in the other candidate nuclei, increasing the chances that a reliable and confirmed result is obtained. Of course, the fact that the nuclear shell model (NSM) consistently predicts M^0ν values that are noticeably smaller than in the other methods remains. Is it the inclusion of complicated configurations (states with high seniority) that causes the reduction of M^0ν in NSM, or is it the inclusion of larger single-particle spaces in most other methods? The answer to this question is obviously needed.

Moreover, there are several more general issues that deserve closer attention of nuclear theorists. One of them is the problem of quenching of the axial weak current matrix elements. That effect is well established in ordinary allowed beta decays, and in the evaluation of the beta strength functions involving low-lying nuclear states. Should an analogous effect be applied to the theoretical evaluation of M^0ν? This is, as yet, an open question. Its solution could affect the true values of M^0ν appreciably, by about 30% according to our estimates.

Problem of quenching is just a part of the determination of effective operators. How does one consistently include the renormalization caused by using only a subset of the general Hilbert space in the calculational framework? Renormalization of the g_A, i.e. inclusion of quenching, is an example of such effect. The problem of taking into account the high momentum part of the nucleon interaction, causing the short range repulsion, belongs to that category as well. There is a consensus now, that the recent developments (see section 6) point us in the right direction. More work in that direction would be clearly beneficial.

Finally, most attention up till now was concentrated on the so-called standard scenario, according to which the 0νββ decay would be caused by the exchange of the light Majorana neutrinos that interact through the canonical left-handed weak currents. Much less attention has been given to the possible mechanism involving heavy, ∼ TeV, particle exchange and thus extremely short range effects. Is it really generally true that in these cases the pionic effects, described in section 5, dominate? If that is the case much more detailed evaluations of the M^0ν in such cases is clearly needed. In particular, in the next few years much progress, one hopes, in the exploration of the TeV mass range particles will be achieved at LHC. If some of the suggested particle physics models (left–right symmetry, R-parity violating supersymmetry) might find support at LHC, more work on the corresponding M^0ν matrix elements will be clearly needed.

The original results reported here were obtained in collaboration with Jonathan Engel, Amand Faessler, Gary Prezeau, Michael Ramsey–Musolf, Vadim Rodin, and Fedor Šimkovic. The fruitful collaboration with them is gratefully acknowledged. The work was supported in part by the US NSF grant 0855538.

The two electron capture decay without neutrino emission requires a special comment. Clearly, when the initial and final states have different energies, the process cannot proceed since energy is not conserved. The radiative process, with bremsstrahlung photon emission, however, can proceed and its rate, unlike all the other neutrinoless processes, increases with decreasing Q value [62]. (However, the estimated decay rates are quite small and lifetimes long). In the extreme case of essentially perfect degeneracy the photon emission is not needed, and a resonance enhancement can occur [63].

In the case of resonance the initial state is the atom (Z, A), stable against ordinary β decay. The final state is the ion (Z − 2, A) with electron vacancies H, H' and, in general, with the nucleus in some excited state of energy E*. The resonance occurs if the final energy $E = E^* + E_H + E_{H^{\prime }}$ is close to the decay Q value, i.e. the difference of the initial and final atomic masses, and a perfect resonance occurs when Q − E is less than the width of the final state which is dominated by the electron hole widths $\Gamma _H, \Gamma _{H^{\prime }}$. The decay rate near resonance is given by the Breit–Wigner type formula

$$\begin{equation} \frac{1}{\tau } = \frac{(\Delta M)^2}{(Q-E)^2 + \Gamma ^2/4}\Gamma , \end{equation} \tag{ A.1 }$$

where ΔM is the matrix element of weak interaction between the two degenerate atomic states.

The states of definite energy, the eigenstates of the total Hamiltonian, are superpositions of the initial and final states, mixed by ΔM. But in reality, the initial state is pure, and not a state of definite energy, since the final state decays essentially immediately.

The mixing matrix element is [63]

$$\begin{equation} \Delta M \sim \frac{G_F^2 \cos ^2 \theta _C}{4 \pi } \langle m_{\beta \beta } \rangle |\psi (0)|^2 g_A^2 M^{0\nu } , \end{equation} \tag{ A.2 }$$

where ψ(0) is the amplitude at the origin of the wave function of the captured electrons and M^0ν is the nuclear matrix element, same one as before. Clearly, if the resonance can be approached, the decay rate would be enhanced by the factor 4/Γ compared to Γ/(E − Q)², where the width Γ is typically tens of eV. Estimates suggest that in such a case the decay lifetime for 〈m_ββ〉 ∼ 1 eV could be of the order of 10^24–25 years. However, chances of finding a case of a perfect (eV size) resonance when E is of order of MeV are not large. Indeed, in the best case found so far, in ¹⁵²Gd [64], the quantity Q − E = 0.91(18) keV, still with the predicted half-life of only ∼10²⁶ years for 〈m_ββ〉 = 1 eV.