The first large-scale shell-model calculation of the two-neutrino double beta decay of $^{76}$Ge to the excited states in $^{76}$Se

Large-scale shell-model calculations were carried out for the half-lives and branching ratios of the $2\nu\beta\beta$ decay of $^{76}$Ge to the ground state and the lowest three excited states $2_1^+$, $0_2^+$ and $2_2^+$ in $^{76}$Se. In total, the wave functions of more than 10,000 intermediate $1^+$ states in $^{76}$As were calculated in a three-step procedure allowing an efficient use of the available computer resources. In the first step, 250 lowest states, below some 5 MeV of excitation energy, were calculated without truncations within a full major shell $0f_{5/2}-1p-0g_{9/2}$ for both protons and neutrons. The wave functions of the rest of the states, up to some 30 MeV, were computed in two more steps by introducing two consecutive stages of truncation. The computed magnitudes of the $2\nu\beta\beta$ nuclear matrix elements (including the value of the axial-vector coupling $g_{\rm A}$), $\vert M_{2\nu}\vert g_{\rm A}^2$, converged to the values 0.168$g_{\rm A}^2$, $1.2\times10^{-3}$$g_{\rm A}^2$, 0.121$g_{\rm A}^2$, and $3.1\times10^{-3}$$g_{\rm A}^2$ for the $0^+_{\rm g.s.}$, $2^+_1$, $0^+_2$, and $2^+_2$ states, respectively. Using up-to-date phase-space integrals, the corresponding branching ratios were derived to be 99.926\%, 4.4$\times10^{-5}$\%, 0.074\% and 2.5$\times10^{-7}$\%. The experimental half-life $(1.926\pm0.094)\times10^{21}$ yr of the ground-state transition was used to derive the value $g_{\rm A}=0.80\pm0.01$ for the axial-vector coupling, which is consistent with other shell-model calculations suggesting a quenched value of $g_{\rm A}$. Using this value of $g_{\rm A}$, predictions for the transition half-lives were derived.

physics, it is interesting in terms of nuclear structure [2], in terms of accessing the effective value of the weak axial coupling g A [3,4], and in terms of being an irreducible background component in 0νββ measurements [5]. In the history of double-beta-decay measurements and calculations 76 Ge, decaying to the ground state and excited states of 76 Se, has been one of the flagship cases. In this Letter, we present the first large-scale shell-model calculation of the nuclear matrix elements (NMEs) of the 2νββ decay of 76 Ge to the 0 + ground state (0 + g.s. ) and 2 + 1 , 0 + 2 and 2 + 2 excited states in 76 Se. Experimentally, the double-beta-decay search for 76 Ge is performed with HPGe semiconductor detectors as suggested by Fiorini et al. [6]. Due to the superb energy resolution of such detectors, it is convenient to use them for searches of peak-like spectra, like is the case for the 0νββ decay of 76 Ge. The two dominant double-beta-decay modes of 76 Ge are the 2νββ decay with the emission of two electrons providing a continuous energy spectrum while the neutrinoless decay mode results in a peak at the Q value. The Q value has been measured to be 2039.061(7) keV [7]. To observe such rare event as double beta decays requires a very low background and various components contribute. Among them are the installation of detectors deep underground, material screening to deduce the background level especially from the natural decay chains of U and Th. Isotopic enrichment in 76 Ge is nowadays mandatory. Further techniques like pulse-shape analysis and veto systems against muons help to reduce background further. The driving experiments in the last decade have been the MAJORANA demonstrator and the GERDA experiment, employing bare HPGe detectors, the latter within liquid argon (see the review [8]). The GERDA has recently measured the 2νββ-decay half-life of 76 Ge to the ground state of 76 Se to be T ββ 1/2 = (1.926 ± 0.094) × 10 21 yr [9]. On the nuclear-structure side, the 2νββ decay of 76 Ge to the ground and excited states has attracted attention mainly in the nuclear-structure community specialized to many-body frameworks based on the quasiparticle random-phase approximation (QRPA) [10]. The ground-state transition has usually been dealt with by the use of the proton-neutron QRPA (pnQRPA) [11] and the excited-state transitions have been handled by the use of various extensions of the QRPA and pnQRPA [12]. In Civitarese et al. [13] foundations for the multiple commutator model (MCM) were laid, and further application of the theory to many 2νββ decays was performed in Aunola et al. [14]. Various higher-QRPA frameworks were applied to 2νββ decays by Bobyk et al. [15], Stoica et al. [16] and Raduta et al. [17]. After the introduction of the renormalised QRPA (RQRPA) by Toivanen et al. [18], the theory was applied to 2νββ decays by the same authors in [19] and by Schwieger et al. in [20]. No shell-model calculations before the present one could be achieved due to limited computational resources which could not match the number of active valence particles in a realistic valence space covering the 0f 5/2 − 1p − 0g 9/2 major shell, for both protons and neutrons. In particular, the involved nuclei, 76 Ge, 76 Se and the intermediate nucleus 76 As, are almost in the middle of the 0f 5/2 − 1p − 0g 9/2 major shell in terms of the proton (Z = 32 − 34) and neutron (N = 42 − 44) numbers.
The half-life of 2νββ decay is given by [2] t (2ν) where G (2ν) is the phase-space integral which captures the kinematics of the decay. For the evaluation of the phase-space integral we use the values (readily available for the decays to 0 + g.s. , 0 + 2 , and 2 + 1 ) and expressions given in [21]. The quantity M 2ν is the 2νββ NME which encodes the nuclear-structure information and can be written for β − β − decay from the initial 0 + ground state where J is the spin of the final state (0 or 2 in the present case), m e is the electron rest mass, E(1 + m ) − M i is the energy difference between the m th intermediate 1 + state and the ground state of the initial nucleus and Q ββ is the energy released in the decay, i.e., the Q value. The power k assumes the value k = 1 for J = 0 and k = 3 for J = 2 (see the derivation and discussion in [22]). The sum runs over the 1 + states in the intermediate nucleus and the exponent of the denominator is one for the decays to the 0 + states and three for the decays to the 2 + states. Because of the higher exponent, the cumulative matrix elements for the decays to the 2 + states are not only much smaller than for the decays to the 0 + states, but they also converge much faster as functions of the energy of the intermediate 1 + states.
The third quantity in eq. (1) is the axial-vector coupling constant g A . Its free-nucleon value determined from the neutron decay half-life is g A = 1.27. However, when calculations for 2νββ decays have been carried out using the nuclear shell model and other microscopic nuclear-structure models, quenched values of g A ≈ 1 (or even less, see the reviews [3,4]) have been consistently needed for reproducing the experimental half-lives (see also the review [5]). In the present study, we predict branching ratios which are independent of the adopted value of g A as long as we make the reasonable assumption that the needed effective value is the same for all the transitions. While the NMEs are not without their own nuclear-structure-related uncertainties (restricted single-particle model spaces, lacking many-body configurations and three-body forces, etc.), the quenching issue is believed to be related to the impulse approximation (see [23]) leading to the omission of the meson-exchange (or two-body) currents, shown to be quite important by Gysbers et al. [24]. The present shell-model calculations exploit an effective nuclear Hamiltonian where the deficiencies pertaining to nuclear structure and lack of meson-exchange currents have been (at least partially) compensated by fitting the two-body Hamiltonian to independent data in a nuclear region relevant for the present calculations. Again, the branching ratios are less sensitive to these deficiencies since they tend to cancel when taking ratios of the transition half-lives.
The nuclear-structure calculations for the wave functions and one-body transition densities (OBTDs), needed for the evaluation of the transition half-lives, were carried out with the shellmodel software NuShellX@MSU [25].
The calculations were carried out in a single-particle model space consisting of the orbitals 0f 5/2 − 1p − 0g 9/2 using the effective Hamiltonian JUN45 [26] . The same Hamiltonian was used to successfully describe observables such as level schemes and half-lives of the neighboring nuclei 71 Ge and 71 Ga in [27]. The computational burden of carrying out the calculations in the entire model space is enormous and even using a powerful computer cluster some compromises had to be made. In any case, we managed to calculate 250 intermediate 1 + states in 76 As using the full model space, as well as the initial state 0 + g.s. in 76 Ge and the four final states of interest 0 + g.s. , 0 + 2 , 2 + 1 and 2 + 2 in 76 Se. There are a couple additional states in 76 Se to which a 2νββ transition is not energetically forbidden but due to the much smaller Q-values and/or larger differences in angular momentum between the initial and final states, the branching ratios to these states are orders of magnitude smaller and can not be experimentally detected in the foreseeable future. The shellmodel Hamiltonian managed to reproduce the low-energy spectrum of 76 Se well: the ordering of the four final states was correct, with the predicted energies of the excited 2 + 1 , 0 + 2 and 2 + 2 states being 562 keV, 1060 keV and 1397 keV, agreeing nicely with the corresponding experimental energies of 559 keV, 1122 keV and 1216 keV [28].
The density of the 1 + states in 76 As is high and we ended up in a three-step shell-model calculation in order to handle these states up to excitation energies that guarantee the convergence of the 2νββ NMEs. In the first step, we managed to calculate the lowest 250 intermediate 1 + states in the full 0f 5/2 − 1p − 0g 9/2 model space, reaching 4.7 MeV in excitation energy. Since we can expect contributions from higher intermediate states to be significant, we carried out additional calculations for these states using two sequential truncations of the model space. With the first truncation, i.e. keeping the neutron orbital 0f 5/2 filled, we managed to calculate also the next 250 1 + states, having now access to the lowest 500 intermediate states in the second step of the calculations. In the third step we kept, in addition, the proton orbital 0g 9/2 empty, which allowed us to calculate the lowest 9999 intermediate states and reach 30.6 MeV in excitation energy. These three steps were combined so that the OBTDs and excitation energies were extracted as follows: For the lowest 250 intermediate 1 + states we used the first-step results of the full-models-space calculation, for the next 250 1 + states we used the results of the second-step truncated calculation and for the rest of the 1 + states we used the results of the third-step truncated calculation. Since the density of states was highest for the first-step full calculation and lowest for the most truncated third-step calculation, we ended up with 10,266 intermediate 1 + states instead of the 9999 states of the most truncated calculation.
Since kinematics have a huge effect on the predicted half-lives, we used the available experimental information to fix the excitation energies of the lowest intermediate 1 + states. The calculations predict 2 − as the ground state in 76 As, in agreement with the data. The computed excitation energies of the lowest three 1 + states are 185 keV, 360 keV and 519 keV, consistently by some 140 − 400 keV higher than the experimental energies of 44.425(1) keV, 87 keV and 120 keV [28]. There is currently no reliable experimental information on the excitation energies of the higher 1 + states. In order to achieve the maximal experimental input, the computed excitation energies of the mentioned 1 + states were shifted for all calculations so as to match the measured energies.
The cumulative nuclear matrix elements M 2ν for the decays to the 0 + g.s. and 0 + 2 states are presented in Fig. 1 and those to the 2 + 1 and 2 + 2 states are given in Fig. 2. The matrix elements for the decays to the 0 + states seem to have converged by 10 MeV in excitation energy and for the decays to the 2 + states by about 3 MeV in excitation energy.
The three-step shell-model application predicts the values of the 2νββ NMEs (including the value of the axial-vector coupling g A ), |M 2ν |g 2 A , to be 0.168g 2 A , 1.2 × 10 −3 g 2 A , 0.121g 2 A , and 3.1 × 10 −3 g 2 A for the 0 + g.s. , 2 + 1 , 0 + 2 and 2 + 2 states, respectively. When only using the first-step results, i.e. the 250 states from the full-model-space calculations, the corresponding matrix elements are 0.166g 2 A , 1.2 × 10 −3 g 2 A , 0.128g 2 A , and 3.1 × 10 −3 g 2 A , respectively. This means that the effects of the higher states on the final NMEs are fairly modest, with the main difference being that the NME for the decay to the 0 + 2 state is somewhat smaller in the three-step calculation, leading to a reduction of the branching ratio to the 0 + 2 state by 12% in relative terms. The second step of the calculations converges the NMEs to essentially their final values as indicated by the dotted lines in Figs. 1 and 2. The final convergence is achieved in step 3. Concerning the accuracy of the present results, one can compare the present results with those of the large-scale jun45 shell-model calculation of Caurier et al. [29] in the same single-particle valence space. The calculation of [29] gives for the ground-state NME the value 0.170g 2 A which is almost exactly the value obtained in the present calculation. Both of these results are also in line with the result of the jun45 calculation of Brown et al. [30]. Furthermore, one has to point out that the presently used single-particle valence space,    [3][4][5] for the decay of 76 Ge to the ground state and lowest three excited states (column 1) of 76 Se. The phase-space integrals for the lowest three states are taken from [21] and the fourth one has been calculated using the formulas of this reference. Experimental Q values (column 2) are used in the calculations. The half-life for the ground-state transition matches the measured one [9] and the rest of the half-lives are based on the shell model calculations using gA = 0.80 ± 0.01, derived from the comparison of the computed half-life with the experimental one for the ground-state transition. The uncertainties include only the uncertainty related to the value of the experimental half-life of the ground-state transition. In the sixth column we list the measured lower limits for the half-lives, including the corresponding reference in the last column. 3.1×10 −3 3.33 × 10 −25 (7.57 ± 0.37) × 10 29 > 1.4 × 10 21 [32] one major shell, is rather limited and contributions from the outside of the valence space may be expected. However, these outside contributions can be, to a large extent, taken into account by the spectroscopy-fitted effective Hamiltonian used in our calculations.
The computed NMEs and phase-space integrals can be used to derive the branching ratios 99.926%, 4.4×10 −5 %, 0.074%, and 2.5×10 −7 % for the transitions to the 0 + g.s. , 2 + 1 , 0 + 2 , and 2 + 2 states, respectively. These branching ratios are shown schematically in Fig. 3. The experimental half-life of (1.926 ± 0.094) × 10 21 yr, obtained for the ground-state-to-ground-state decay [9], can be achieved with the value g A = 0.80 ± 0.01 of the axial coupling, which is in line with what has been seen in many shell-model calculations, i.e., that a quenched value of the axial coupling is needed to reproduce the experimental results (see the review [3]). Using this value of the axial coupling, the half-lives for the excited-state transitions listed in Table 1 were derived. The uncertainties of the computed half-lives include only the uncertainty related to the value of g A , i.e., to the experimental half-life of the ground-state transition. This is because the uncertainties related to the matrix elements cannot be reliably estimated, though our adopted method controls for systematic underor overestimation related to the valence space by fixing the ground-state-to-ground-state half-life using experimental data. In the same table the currently known experimental half-lives and their lower limits are shown. As can be seen, the experimental lower limit for the half-life of the transition to the 0 + 2 state is some two orders of magnitude lower than the computed half-life. For the other transitions the experimental half-life limits are way lower than the computed half-lives.
In Table 2 we compare the presently obtained shell-model results for the ground-state and excited-state transitions with the corresponding results of other calculations found in the literature (excluding the already mentioned two shell-model results [29,30] for the ground-state transition). For the ground-state transition the presently computed NME, 0.168, reproduces the measured half-life with g A = 0.80 ± 0.01 and the rest of the calculations with quenched values in the range g A = 0.97 − 1.2. For the decay to the 2 + 1 state the presently obtained result is quite close to the MCM result of [13] and the HRPA result of [17]. Also the other results are not very far from the present result, except the HQRPA result of [16]. For the decay to the 0 + 2 state the presently computed NME is quite close to the RQRPA one of [19] and not very far from the MCM result of [13] and the HQRPA result of [16]. Concerning the decay to the 2 + 2 state, the presently obtained NME is consistent with the MCM results of [13,14]. Overall, the present results for the excited-state NMEs correspond surprisingly closely to the MCM results of Civitarese et al. [13].  In this Letter we present the first large-scale shell-model calculations concerning the 2νββ decay of 76 Ge to the ground state and the lowest three excited states 2 + 1 , 0 + 2 and 2 + 2 in 76 Se. Due to the fact that the proton and neutron Fermi surfaces of the involved nuclei are roughly at the middle of the 0f 5/2 − 1p − 0g 9/2 major shell, adopted as the valence space of our calculations, the computational challenges are formidable even for a powerful computer cluster used in the present calculations. This is why we performed the computations of the related nuclear matrix elements by introducing a three-step procedure. The first step allowed us to extract the wave functions of the intermediate 1 + states in 76 As in a full valence space below the low-energy limit of some 5 MeV. These states essentially determine the magnitude of the NMEs for the 2 + states and carry the bulk of the contributions to the NMEs of the 0 + states. The wave functions for the states above the low-energy limit were computed in two steps by imposing two consecutive truncations in the valence space. All the NMEs were converged very close to their final values already after the first truncation in step 2. The third step guaranteed the convergence of the NMEs, and they were subsequently used to derive branching ratios for the decays to the ground and excited states by adopting up-to-date phase-space integrals. Concerning decay transitions to these excited states, our computed branching ratios can be used as guidelines in designing future 2νββ-decay experiments using a 76 Ge source.
The measured half-life of the ground-state transition could be achieved by using the value g A = 0.80 ± 0.01 for the axial-vector coupling. Adopting this value of the axial coupling, the halflives for the decay transitions to the ground and excited states were derived. Comparing our results for the NMEs with those of other calculations, based on higher QRPA schemes of various kinds, showed a reasonable overall agreement with most of the calculations for all NMEs. In particular, we recorded a very good agreement with the MCM results of Civitarese et al. [13] for the excited-state NMEs. The step-wise shell-model calculations, introduced by us in this Letter, could possibly be used in other computationally intense 2νββ-decay and/or 0νββ-decay calculations in the future. An other possible way to speed up the calculations and possibly avoid model-space truncations is to use the Lanczos algorithm as described in [33,34,35].