Neutrinoless ββ decays to excited 0 + states and the Majorana-neutrino mass

The nuclear matrix elements (NMEs) corresponding to the neutrinoless double-β (0 νββ ) decays to excited 0 + states of major experimental interest are calculated. All these decay transitions are electron emitting (0 νβ − β − decays) and take place in the mass A = 76 , 82 , 96 , 100 , 110 , 116 , 124 , 130 , 136 nuclei. This work is an extension of our previous work [Phys. Rev. C 91 , 024613 (2015)], where 0 νββ decays to the ground states of the same nuclei were treated. We calculate the NMEs for transitions mediated by both the light (l-NMEs) and the heavy (h-NMEs) Majorana neutrinos. A higher-QRPA (quasiparticle random-phase approximation) framework, the multiple-commutator model, is adopted for the calculations, including a previously omitted contribution to the transitions to two-phonon states. A Bonn G -matrix-based effective nucleon-nucleon interaction is generated by exploiting the recently proposed isoscalar-isovector decomposition of the particle-particle proton-neutron interaction parameter, g pp . All the appropriate short-range correlations, nucleon form factors, and higher-order nucleonic weak currents are included to benchmark our calculations. The relevant nuclear spectroscopy was checked to validate the nuclear models used. The computed l-NMEs and h-NMEs are compared with the available other calculations and the relevance of the new included two-phonon term is discussed. The results are summarized by easy-to-use half-life-Majorana-mass interrelations.


I. INTRODUCTION
The neutrino-oscillation experiments have succeeded in gathering precise information on neutrino-mass differences and neutrino-mixing amplitudes during the past decade [1]. The neutrinoless double-β (0νββ) decay serves as a complementary means of extracting information on the fundamental properties of neutrino, like its basic character as a Dirac or Majorana particle [2][3][4][5]. Because the 0νββ decays take place in atomic nuclei one needs the aid of nuclear-structure physics to extract quantitative information about neutrino properties from the results of the 0νββ experiments. This is where the calculation of the associated nuclear matrix elements (NMEs) comes into play. Accurate calculation of the NMEs is a challenging task and a central issue in the present-day double-β-decay community. A good selection of the presently used nuclear models and their 0νββ results can be found in the review [6].
Until quite recently most of the 0νββ efforts have concentrated on the electron-emitting (0νβ − β − ) decays to the ground state, 0 + gs , of the daughter nucleus [2]. In the case of the two-neutrino double-β (2νββ) decay the first experimental results on the transitions to the first excited 0 + state, 0 + 1 , in 100 Mo and 150 Nd have been obtained (see the recent review [7]). On the theory side, first attempts to compute the 2νββ-decay rates to excited states in the quasiparticle random-phase approximation (QRPA) formalism were made in Refs. [8][9][10][11]. Recent theoretical results and references to earlier ones can be found in Ref. [12].
Until now, no positive signal has been detected in the 0νββ experiments, in particular concerning decays to excited final states. Attempts have been made to compute the NMEs related to 0νββ decays to excited 0 + states by using the QRPA framework [13][14][15][16][17][18][19], the interacting shell model (ISM) [20], and the interacting boson approximation (IBA-2) [21]. In the QRPA case the basic QRPA was extended to the multiple-commutator model (MCM) [10,22] by combining the charge-conserving QRPA (ccQRPA, for the excited states in an even-even nucleus) with the proton-neutron QRPA (pnQRPA, for the states of an odd-odd nucleus) using a common even-even reference ground state. Decay transitions to the one-and two-ccQRPA-phonon states were considered and the corresponding transition amplitudes were written explicitly in Ref. [10]. Earlier, a slightly different formalism, the two-phonon recoupling approach [8], was used in the description of the 100 Mo decay to the two-phonon 0 + state in 100 Ru.
In the present work we extend the previous MCM calculations by computing 0νβ − β − decays to excited one-and two-phonon states for all nuclei of major experimental interest, including also the heavy Majorana-neutrino exchange as a possible mechanism to mediate 0νββ transitions. At the same time we extend the two-phonon formalism of Ref. [10] by adding an extra term in the decay amplitude, previously omitted as negligible for the light-Majorana-neutrino exchange but possibly important for the exchange of a heavy Majorana neutrino. This work is a continuation of our previous work [23] on 0νββ transitions to the ground state by using the efficient recursive numerical methods reported in Ref. [24]. Like in Ref. [23], we adopt the isospin-symmetry restoration scheme introduced in Ref. [25]. In this scheme one divides the particle-particle proton-neutron interaction constant, g pp [26], into isoscalar (T = 0) and isovector (T = 1) parts and fits the isovector part such that the magnitude of the Fermi NME of 2νββ decay is pushed to zero, thus restoring the isospin symmetry. The isoscalar part of g pp we treat as discussed in Ref. [23]. Like in Ref. [23], we take into account the appropriate short-range nucleon-nucleon correlations [27][28][29] and contributions arising from the induced currents and the finite nucleon size [30].
This article is organized as follows. In Sec. II we introduce the MCM expressions for the one-and two-phonon transition amplitudes. In Sec. III we display and discuss the obtained double-β-decay results for the l-NMEs and h-NMEs. The final conclusions are drawn in Sec. IV.

II. MCM TRANSITION AMPLITUDES
The computational formalism to obtain the 0νββ-decay amplitudes for the transitions to the ground state was thoroughly discussed in Refs. [23,24] and we refer the reader to these sources for further details. Here we display only those expressions that are essential for understanding the background of the present calculations of transitions to the excited 0 + states.

A. Structure of the 0νββ-decay NME
We start by introducing the 0νββ-decay NME that was written in Eq. (38) of Ref. [23] for the ground-state-to-groundstate decays. Its extension to ground-state-to-excited-state decays can be written in the form where k 1 and k 2 label the different pnQRPA solutions for a given multipole J π , stemming from the parent and daughter nuclei of the 0νββ decay. The operators O K in the reduced twoparticle matrix element denote the Fermi (K = F), Gamow-Teller (K = GT), and tensor (K = T) parts of the double-β operator, given explicitly in Ref. [23]. The two-particle matrix element contains also the appropriate short-range correlations, higher-order nucleonic weak currents, and nucleon form factors [23]. The last line of (1) contains the one-body transition density between the initial ground state (0 + i ) and the intermediate state J π k 2 , and it can be obtained from a pnQRPA calculation as discussed in Sec. II C of Ref. [23]. The one-body transition density between the intermediate state J π k 2 and the final 0 + state (0 + f ) can be evaluated for the 0 + ground state by the procedure outlined in Sec. II C of Ref. [23]. To evaluate it for the transitions to the excited 0 + states, one has to resort to the higher-QRPA MCM formalism presented in the following section. The term between the one-body transition densities is the overlap between the two sets of intermediate states emerging from the two pnQRPA calculations based on the parent and daughter even-even ground states.
A possible serious shortcoming of the present QRPA calculations, and also of many other nuclear-structure calculations in different model frameworks, is that they may contain spurious center-of-mass (c.m.) contributions which contaminate the computed nuclear wave functions. Because the c.m. momentum operator is a vector, of spin-parity J π = 1 − , the spurious components tend to affect the computed 1 − states. For general considerations of restoring the Galilean invariance in nuclear many-body problem, see Ref. [31] and the references therein. There exists an exact method for the RPA calculations to remove the c.m. contributions from the wave functions [32] and this method can be probably generalized to the QRPA case. However, we leave a more quantitative analysis considering these c.m. effects for the future.

B. Decay amplitudes in the MCM formalism
Let us write the solution of the pnQRPA equations (see, e.g., Ref. [33]), where we use the shorthand ω = J π k for the kth intermediate state of spin parity J π . Here |QRPA is the QRPA ground state. However, the first excited 2 + state, 2 + 1 , in an even-even (daughter) nucleus is described in the ccQRPA formalism, and the corresponding wave function can be presented as (see [33]) where the normalized two-quasiparticle operators are defined as for any state of angular momentum I in the even-even nucleus. We denote alsoÃ ab (I M) ≡ (−1) I +M A ab (I,−M). It should be noted that here the summation over a b guarantees that there is no double counting of two-quasiparticle configurations and this with the normalized operators (4) guarantees that the wave function is properly normalized with the normalization condition [33] a b The first 2 + state of an even-even vibrational nucleus (which we assume our 0νββ daughter nuclei are) is usually very collective [33] and we refer to the creation operator Q † (2 + 1 ,M) of (3) as a creation operator for a ccQRPA phonon.
For calculational convenience it is preferable to go from the restricted sum of (3) to a nonrestricted one by introducing the correspondence where the barred two-quasiparticle operators are the ones of (4) without the normalizer N ab (I = 2). Then the normalization 064306-2 condition of the 2 + 1 state reads At the same time the two kinds of X and Y amplitudes are related bȳ for any ω = I π n , n denoting the nth even-even daughter state of spin-parity I π . The barred amplitudes are symmetrized ones and possess the convenient symmetry relations (to generate amplitudes with a > b) (10) Using the above-defined barred operators, one can write the two-phonon states of the even-even daughter nucleus in the form Let us denote the one-phonon states in the 0νββ daughter nucleus in general by with ω f = I π n identifying the nth I π state. Then, in the MCM framework the one-body transition densities from the intermediate |ω M state (2) to the final one-phonon state (12) and two-phonon state (11) are calculated by first writing the transition densities as ground-state-averaged multiple commutators and then applying the quasiboson approximation [33] by replacing the QRPA vacuum by the BCS vacuum when taking the ground-state average. The averaged multiple commutators then become where | is the BCS ground state and we have denoted the β − type of transition density operator by with c † p creating a proton on orbital p andc n annihilating a neutron on orbital n.
To convert (13) and (14) to reduced transition densities, needed in Eq. (1), we can use the Wigner-Eckart theorem [33]. We then obtain after a lengthy but straightforward calculation the needed one-body transition densities for the one-phonon final state and for the two-phonon final state with The auxiliary geometric factors read K (2) I LJ = (−1) J +L+j n +j n +1 × K (1) I LJ (j p → j n ; j p → j n ; j n ↔ j p ) = (−1) J +I +j p +j n 2 2 I j n j n j n J L I j n j n j p .
In the earlier works [10,[13][14][15][16][17][18][19]34] the latter term, T (2)ω I LJ , of the two-phonon transition density (18) was omitted owing to its potentially negligible contribution to the β-decay and 064306-3 double-β-decay rates. Here we take full account of (18) to study quantitatively the contribution of the omitted term to the 0νββ decays mediated by a light or a heavy Majorana neutrino.

A. Model parameters and 0 + excited states
The model parameters were fixed by the procedures advanced in detail in our previous work [23], in Sec. III A. There the construction of the various single-particle model spaces with their single-particle energies is thoroughly explained. For the pnQRPA the particle-hole parameter g ph and the particle-particle parameter g pp are adopted. The g pp parameter is further decomposed into isoscalar (T = 0) and isovector (T = 1) parts and the isovector part is fixed such that the magnitude of the Fermi NME of 2νββ decay becomes zero. This guarantees the restoration of the isospin symmetry [25]. The isoscalar part of g pp was fixed in a way discussed in Ref. [23].
To describe the excited states of the 0νββ daughter nuclei, we use the ccQRPA with the associated particle-hole parameter G ph , and the particle-particle parameter G pp . For the 2 + mode the G ph parameter we fixed by reproducing the excitation energy of the first (collective) 2 + state in the ccQRPA calculation. For the particle-particle parameter we used the default value G pp = 1.0 because the value of this parameter has practically no effect on the properties of the 2 + 1 state used in the two-phonon calculations. The 0 + state, used in the calculations of decay transitions to the one-phonon states, was fitted to its experimental energy by varying simultaneously the G ph and G pp parameters, removing at the same time the spurious lowest 0 + root of the ccQRPA by forcing its energy to zero [35]. The same nucleon-nucleon interaction (Bonn G matrix) and BCS vacuum was used in the pnQRPA and the ccQRPA calculations.
In the following we review our results concerning the l-NMEs and h-NMEs of the mass A = 76,82,96,100, 110,116,124,130,136 nuclear systems. In the mass A = 76,82,136 systems the two-phonon 0 + state in the daughter nuclei 76 Se , 82 Kr, and 136 Ba is clearly identifiable [13,17] at energies E(0 + 2-ph ) = 1.1223,1.4876,1.5790 MeV, respectively. In Ref. [16] the mass A = 124,130,136 systems were discussed and the two-phonon states in the daughter nuclei 124 Te and 130 Xe were identified at energies E(0 + 2-ph ) = 1.6573,1.7935 MeV, respectively. In Ref. [14] the two-phonon states in the daughter nuclei 96 Mo and 116 Sn were identified as the second excited 0 + states at energies E(0 + 2-ph ) = 1.330,2.0275 MeV, respectively. At the same time the first excited 0 + states in these nuclei were assumed to be one-phonon vibrational states at energies E(0 + 1-ph ) = 1.1481,1.7569 MeV, respectively. In Refs. [11,15] the first excited 0 + state in 100 Ru was treated as a two-phonon state [E(0 + 2-ph ) = 1.1303 MeV] and the second excited 0 + state was treated as a one-phonon vibrational state [E(0 + 1-ph ) = 1.7410 MeV]. Finally, the first excited 0 + state in 110 Cd is a very likely two-phonon excitation at an energy E(0 + 2-ph ) = 1.4731 MeV (see Ref. [18]). In the present work we use all the aforementioned assignments of energies and characters of the excited 0 + states in the 0νββ daughter nuclei in our computations.

B. Comparison of the calculated nuclear structure with experiment
The purpose of this section is to test the validity of the adopted nuclear-structure models by comparing the calculated quantities with available experimental data. These quantities include the low-energy 0 + and 2 + spectra, the two-phonon energies, and the B(E2) values.

Low-energy 0 + and 2 + spectra
Let us first look at the energy spectra of the low-lying 0 + and 2 + excitations. The two-phonon states, used to model the structure of the low-energy 0 + excitations in this paper, are built from two QRPA quadrupole phonons. Moreover, the first 0 + excitation in nuclei 96 Mo and 116 Sn and the second 0 + excitation in the nucleus 100 Ru are treated as one-phonon states. Thus, comparing the theoretical 0 + and 2 + excitation energies with the corresponding experimental spectra may give some preliminary hint about the validity of the one-and twophonon wave functions used in the computation of the NMEs. Figure 1 presents the theoretical and experimental lowenergy 0 + spectra for the nuclei 96 Mo and 116 Sn. The energy of the 0 + 2 two-phonon state in 96 Mo compares well with the experimental energy of 1.330 MeV, although the experimental value is not certain, as indicated by the question mark (?). The value shown here is from the latest Nuclear Data Sheets [36] analyzing this mass region. In the nucleus 116 Sn the predicted two-phonon state is found at the energy of 2.584 MeV. This overestimates by a half of a MeV the corresponding experimental energy at 2.027 MeV. However, the rest of the 0 + spectrum compares very well with the experimental levels. It seems that low-energy 0 + spectra in 96 Mo and 116 Sn are quite well modeled as one-and two-phonon excitations. Figure 2 shows the theoretical and experimental low-energy 2 + spectra for the nuclei 96 Figure 3 shows the theoretical and experimental low-energy 2 + spectra for the nuclei 82 Kr and 124 Te. In the nucleus 82 Kr the predicted two-phonon state at the energy of 1.550 MeV is in good agreement with the corresponding experimental level at 1.475 MeV. One can identify the higher experimental levels as one-phonon excitations as well, although the multipolarities of the higher levels are not certain. The predicted two-phonon state of the nucleus 124 Te at the energy of 1.206 MeV is also in good agreement with the experimental level at 1.326 MeV. The higher energy experimental spectrum is now much more dense and not given by the one-phonon excitations alone. Figure 4 displays theoretical and experimental low-energy 0 + and 2 + spectra for the nuclei 82 Kr and 130 Xe, respectively. In 82 Kr the two-phonon model predicts a correct energy for the first 0 + excitation. The second excited state at the energy of 2.139 MeV compares favorably with experiment as well. The low-energy 0 + spectrum in 82 Kr seems to have a clear one-and  two-phonon structure, as was also the case for the 2 + spectrum. The two-phonon model also predicts a correct energy for the second 2 + excitation in the nucleus 130 Xe. Identification of the higher experimental states as one-phonon states is not so clear because of the uncertainties in the experimental multipolarities and because of the (thus resulting) higher level density. The experimental level at the energy of 3.230 MeV is given a relatively large error estimate ±0.200 MeV so it might compare favorably with the QRPA state at 3.542 MeV. Figure 5 presents the theoretical and experimental lowenergy 0 + and 2 + spectra for the nucleus 136 Ba. The twophonon model works very well in predicting both the 0 + 1 and the 2 + 2 energies. In case of the higher-energy states, one can identify the QRPA-predicted 2 + one-phonon states at the energies 2.167 and 2.235 MeV with the experimental states at the energies 2.129 and 2.223 MeV, respectively. Of course, one would also need to investigate the electromagnetic properties of these states to confirm this identification. The higher QRPA 0 + spectrum compares nicely with the experimental levels as well. Based on this comparison it seems that the low-energy 0 + and 2 + excitations in 136 Ba are given by one-and two-phonon excitations alone.

Two-phonon energies and B(E2) values
In Table I Table I. The theoretical values in Table I are given by the following very simple systematics, valid for the two-quadrupole-phonon excitations [33],  [45]. Then e (n) eff was adjusted in such a way that the experimental B(E2; 2 + 1 → 0 + gs ) value is reproduced by the theory. The thus-obtained values for the effective charges are listed in the last row of Table I. The rest of the theoretical B(E2) values, listed in Table I in Weisskopf units, then follow from the two-phonon formulas (23) and (24). In Table I we have also compared the theoretical energies and B(E2) values with available experimental data. It can be seen that the experimental energies agree quite well with the energies of the degenerate two-phonon triplet given by Eq. (25). However, in case of the B(E2) values it seems that only the 4 + 1 states follow closely the two-phonon prediction (23). For the states 0 + and 2 + the two-phonon estimate (23) is usually larger than the experimental value by a factor of two. Still, it is remarkable that the experimental B(E2; 2 + 2 → 0 + gs ) values are close to zero and confirm the two-phonon picture where these B(E2) values are exactly zero. Contrary to the general trends the two-phonon picture fails completely in the description of the 0 + 2 → 2 + 1 E2 transition in the nucleus 116 Sn. This could have consequences for the reliability of the 0νββ NMEs discussed in the next sections. One should also mention that the use of the microscopic two-phonon description violates slightly the Pauli principle because both 2 + phonons contain the same intrinsic structures. However, this violation is not expected to affect significantly the final results in two-phonon calculations [46]. Overall, it seems that the two-phonon model of Eqs. (23)- (25) gives pretty good estimates for most of the energies and B(E2) values of relevance in this work.

C. Matrix elements for light-neutrino exchange
We start the presentation of our 0νββ results with a discussion of the NMEs corresponding to the exchange of a light Majorana neutrino (the l-NMEs). The l-NMEs 064306- 6  TABLE II. Values of the computed l-NMEs for 0 + gs −→ 0 + 1 transitions. Columns 3-5 show the decomposition of the presently calculated total NMEs (column 6) in terms of the Fermi, Gamow-Teller, and tensor contributions. Our present total NMEs are compared with the IBM-2 NMEs of Ref. [21] in the second-to-last column and with the earlier QRPA computations of Refs. [14,15,17,18] in the last column.  [17] corresponding to the 0 + gs −→ 0 + 1 decay transitions have been decomposed according to the Fermi, Gamow-Teller, and tensor contributions in Table II. The tensor contribution seems to become negligible for the transitions to the excited states. It has some influence to the ground-state-to-ground-state transitions, although the effect is quite small in that case also [23]. We have compared our present results with the IBM-2 NMEs of Ref. [21] in the second-to-last column of Table II. Agreement between the QRPA and IBM-2 calculations is satisfactory for the lighter nuclei, larger deviations occurring for the heavier nuclei. There seems to be a notable discrepancy in the case of the decay of 96 Zr. The IBM-2 results are computed by using the Miller-Spencer (M-S) parametrization of the Jastrow function which imitates the effect of the hard-core nucleon-nucleon repulsive force on the two-body wave functions. In this work we have used the softer CD-Bonn parametrization. The difference between the NMEs calculated with different short-range correlations is, however, quite small for the decays to the two-phonon excited states [17]. This is one of the different features between the ground-state and the excited-state transitions. For the ground-state transitions the difference between the M-S and CD-Bonn parametrizations is large (of the order of 1.0), whereas for the two-phonon excited-state transitions it is rather small (of the order of 0.1). This property was noticed in Ref. [17], where the effects of the Unitary Correlation Operator Method and Jastrow short-range correlations on the two-phonon transitions were compared. In case of the nucleus 96 Zr there exists a strong SRC dependence similar to the ground-state decays, but this dependence is not even nearly strong enough to suppress the NME to the level of the IBM-2 result.
We have also compared our present results with the earlier QRPA calculations of Refs [14,15,17,18] in the last column of Table II. In these earlier computations the second term (20) of the two-phonon transition density was omitted as negligible. In the calculations of Refs. [14,15] the computational scheme was based on the relativistic harmonic confinement model (RHCM; see Ref. [47]) instead of the more phenomenological standard formulation used in this work. In RHCM, the nucleon current emerges from the quark dynamics inside the nucleons. No short-range correlations are taken into account beyond the RHCM-predicted nucleon form factors. In addition, the coupling parameters g V and g A do not explicitly appear because they are set to unity on the quark level and calculated from the quark dynamics on the nucleon level [2,48]. This is why the comparison of the present results, obtained in the standard formulation, and the results of Refs. [14,15] is not straightforward. It was noticed in Ref. [17] that in the RHCM formulation the decays to the two-phonon states are suppressed more than in the standard formulation. It seems, according to our present calculations, that also the decays to the one-phonon states (decays of 96 Zr and 116 Cd; see Table II) are suppressed more in the RHCM formulation.
In Figs. 6(a)-6(i) we have plotted the QRPA multipole decomposition of the total l-NME M (0ν) for all nuclei considered in this paper and for 0 + gs −→ 0 + 1 decay transitions. All decompositions are calculated with the bare value of the axial-vector coupling g A = 1.26. Matrix elements depend mostly on the details of the 1 + multipole. In some nuclei also the 0 + , 2 + , and 2 − contributions play a significant role. Structure of the multipole distribution seems to be the same for one-and two-phonon transitions.
The l-NMEs corresponding to the 0 + gs −→ 0 + 2 transitions have been decomposed according to the Fermi, Gamow-Teller, and tensor contributions in Table III. The l-NMEs for these 064306-7 transitions are usually suppressed by one order of magnitude relative to the 0 + gs −→ 0 + 1 matrix elements and because of this are probably not so interesting from the point of view of experimental detection of the 0νββ.
In Fig. 7 we have analyzed the effect of the previously omitted second term (20) of the two-phonon transition density on the total l-NMEs. For the heavier nuclei ( 116 Cd , 124 Sn , 130 Te, and 136 Xe) the shift caused by the second term is quite small 064306-8  (20) included. It is seen from the figures that the new term basically switches on the J = 0 contribution, which becomes comparable in magnitude with the earlier dominating J = 2 contribution. The two angular-momentum contributions interfere destructively, which causes the suppression of the matrix elements. Similar cancellations between the J = 0 and J = 2 contributions are observed also in the shell-model calculations [49], although in the shell model the effect seems to be more intense, leading to more suppressed NMEs. In Figs. 9(a) and 9(b) we have plotted the pair decompositions also for the first-excited-state transitions in nuclei 96 Zr and 116 Cd. In these transitions the final state is a one-phonon vibrational state as described in Sec. III A. The structure of the pair decomposition is more simple in the one-phonon case than in the two-phonon case, the J = 0 part clearly dominating over the rest of the contributions. The destructive interference effect is also absent, leading to potentially larger matrix elements.
The 0νββ transition operators O K in the matrix element expression (1) depend the relative radial coordinate r between the two decaying neutrons (see the explicit expressions for the transition operators in Ref. [23]). In Figs. 10(a) and 10(b) we have plotted the radial dependence of the matrix elements M (0ν) F , M (0ν) GT , and the full NME M (0ν) for the nuclei 76 Ge and 96 Zr. The tensor matrix element was omitted owing to its negligible size. The curves in the figures correspond to transitions to the first excited 0 + states. Complete matrix elements are found by integrating over the radial curves The two sets of curves shown have very different shapes. This difference can be traced back to the different structure of the decay final states. For 76 Ge the final state has a two-phonon structure, whereas for 96 Zr the final state is a basic one-phonon excitation. The radial curves for 76 Ge reveal why the two-phonon transitions do not show strong dependence on the parametrization of the short-range-correlation effects. A large contribution to the NME comes from higher relative distances, so the SRC effects play only a minor role. By contrast, for the ground-state-to-ground-state transitions most of the NME seems to come from shorter relative distances r 2 fm [23,50], which leads to a strong SRC dependence for those decays. The radial curves for 96 Zr resemble more those found for the ground-state decays. Most of the matrix element comes from short relative distances leading again to a more intense SRC dependence.

D. Matrix elements for heavy-neutrino exchange
Next we present our results for the NMEs corresponding to the exchange of a heavy Majorana neutrino (the h-NMEs). The h-NMEs corresponding to the 0 + gs −→ 0 + 1 decay transitions are decomposed according to the Fermi, Gamow-Teller, and tensor contributions in Table IV. Unlike the l-NMEs, the h-NMEs depend strongly on the parametrization of the short-range correlation effects. In this respect they behave similarly to the ground-state transitions. Because of this, we have computed the h-NMEs also with the M-S parametrization to facilitate a more straightforward comparison with the M-S parametrized IBM-2 NMEs of Ref. [21], displayed in the last column of Table IV. In addition to this, we have also included in Table IV the h-NMEs computed without the second term (20) of the two-phonon transition density and by using the M-S SRC. Present QRPA results seem to be only in a partial agreement with the IBM-2 calculations. There exists notable differences for the nuclei 96 Zr and 110 Pd. It was observed in our previous work [23] that for the ground-state decays the IBM-2 h-NMEs are always, approximately, by a factor of 2 smaller than the corresponding QRPA h-NMEs. This systematic deviation seems to be lost for transitions to the excited states. In some cases the IBM-2 NMEs are larger than the corresponding QRPA NMEs in other cases they are smaller, as one notices from Table IV. 064306-9    (20) of the two-phonon transition density is omitted. Our present results are compared with the IBM-2 calculations [21] in the last column.

Nuclear transition
g A h-NMEs, present results M (0ν) [21] M (0ν) (20) We observe that the h-NME for 110 Pd suppresses dramatically when the SRC parametrization is changed from CD-Bonn to M-S. For the majority of nuclei the change from CD-Bonn to M-S SRC cuts the matrix element about 50%, which is a similar suppression as in the case of the ground-state transitions [23]. The effect of the added second term (20) of the two-phonon transition density seems to be much stronger in the heavy-neutrino exchange than it was in the light-neutrino exchange. We have investigated this effect in Fig. 11. One notices from the figure that the second term gives very large shifts to the matrix elements of the lighter nuclei and its omission is no longer acceptable. In fact, for the nuclei 100 Mo and 110 Pd the negative contribution given by the term (20) is so large that it changes the signs of the total NMEs for these nuclei. In Table IV we have changed their sign again to make all the total NMEs positive. The overall sign of the NMEs does not matter if the S-wave mass mode is considered as the only possible decay mechanism as is done in this paper. The sign has to be taken into account, however, if one wishes to compute the matrix elements and half-life contributions coming from the right-handed currents and electron P waves. The nuclei in the decay chain of 136 Xe behave exceptionally because adding the second term increases the magnitude of the NME, whereas for the rest of the nuclei the NME is suppressed. The effect of the added term (20) is, as in the case of the light-neutrino exchange, most easily grasped by observing the pair-angular-momentum decompositions. In Figs  nuclei 76 Ge , 82 Se , 110 Pd, and 136 Xe. As in the light neutrino exchange, the added term (20) switches on the J = 0 contribution, but now the amplification is much larger than in the previous case. In the nucleus 76 Ge the J = 0 term almost exactly cancels the earlier dominating J = 2 contribution. The same thing happens in the nucleus 110 Pd. In 110 Pd, the destructive interference of the J = 0 and J = 2 terms is even more intensified if M-S parametrization is used for the SRC effects, as one notices from the sixth and seventh columns of Table IV. In the nucleus 136 Xe the amplification of the negative J = 0 contribution is milder, which causes the +38% increase of the NME observed earlier.
In Figs. 13(a)-13(d) and 14(a)-14(e) we have plotted the QRPA multipole decomposition of the h-NME for all the nuclei considered in this paper and for the 0 + gs −→ 0 + 1 decay transitions. Multipole decompositions have more irregular shape for the heavy-neutrino exchange than for the light-neutrino exchange. Multipole strength is distributed towards higher angular-momentum values. There exists no characteristic leading multipole either, as was 1 + in the light-neutrino exchange. The multipole structure of 96 Zr differs somewhat from the others. It is more reminiscent of the ground-state-to-ground-state multipole distributions [23].
In Table V we have presented our results for the h-NMEs for the 0 gs −→ 0 + 2 decay transitions. For the nucleus 96 Zr the transition to the second excited state is suppressed by two orders of magnitude relative to the transition to the first excited state, whereas an order-of-magnitude suppression occurs for the nuclei 116 Cd and 100 Mo. The tensor contribution seems to become very important for the decay of 100 Mo.

E. Half-lives for the light-and heavy-neutrino exchange
Starting from the expression (13) of Ref. [23] the 0νββ half-life for the light-neutrino-exchange mechanism can be written as where the effective neutrino mass should be given in units of eV. Similarly, we find for the half-life of the decay transitions mediated by heavy electron neutrinos, the expression where the effective inverse mass of the heavy neutrino, η M , has been defined, e.g., in Ref. [4] [Eq. (66)] or in Ref. [23] and it should be inserted in units of inverse GeV. Although we have not considered the possible particle physics mechanisms for the generation of the neutrino masses in this paper, concentrating on the nuclear-structure aspects of the problem instead, it is worth pointing out that in Ref. [51] the heavy neutrino contribution was shown to be negligible in the context of the seesaw type I scenario. It is likewise appropriate to mention that the h-NMEs contribute also to the processes triggered by the right-handed currents (see Eq. (67) of Ref. [4]).
We have listed the nuclear-structure coefficients C 0ν l and C 0ν h in Table VI. To contrast with the corresponding coefficients for the ground-state-to-ground-state transitions, we have also included the ground-state coefficients as extracted from our earlier results in Ref. [23]. The ranges of the coefficients C 0ν l and C 0ν h correspond to the range g A = 1.00-1.26 for the axial-vector coupling constant. It is seen that for the light-neutrino exchange the half-lives for the transitions to the first excited states are usually an order of magnitude longer than for the ground-state-to-ground-state transitions. In some cases the difference is two or even three orders of magnitude, like for the nucleus 110 Pd. The large suppression of the decay   of 110 Pd is related to the combined effect of a small NME and an unfavorable phase space for this decay. The decay of 96 Zr to the first excited 0 + state in 96 Mo seems to be enhanced relative to the other excited-state transitions. One notices from Table VI that this transition in 96 Zr is not very much suppressed relative to the ground-state transition. This effect was already noticed in Ref. [14]. It was also discussed in Ref. [14] how this enhancement could be utilized in the experimental search for the 0νββ decay. The main benefit is the possible reduction of the background of the measured spectra by the coincidence of the two emitted electrons with the signals emerging from the γ decay of the excited final state. Although 96 Zr seems to be a very promising candidate for these kind of experimental investigations, one must remember that our present QRPA results for 96 Zr are not supported by the IBM-2 calculations, which give an exceptionally small matrix element for this decay (see Table II). Further studies of the excitedstate transitions are needed to settle this discrepancy. Another good option for the experimental detection of an excited-state transition would be 136 Xe, at least from a theoretical point of view.
For the heavy-neutrino exchange the half-lives of the excited-state transitions are usually much longer than for the ground-state transitions. This is particularly striking for the transition to the 0 + 1 state in 110 Cd and transitions to the 0 + 2 states in 96 Mo and 116 Sn. Contrary to this, only an order-of-magnitude suppression occurs for the transition to the 0 + 1 state in 96 Mo. This is the strong 0 + vibrational state enhancing also the corresponding transition mediated by the light Majorana neutrino.

IV. CONCLUSIONS
As an extension to our previous work on the ground-stateto-ground-state transitions, in this work we have calculated the nuclear matrix elements of the neutrinoless double-β-minus decays, mediated by the light or heavy Majorana neutrino, to excited 0 + final states. The matrix elements have been computed for those decay transitions which potentially attract current and future experimental interest. We have decomposed our matrix elements in different ways to allow transparent comparison with other calculations using different nuclearstructure frameworks.
Our calculations have been done by using realistic two-body interactions and single-particle bases. All the appropriate short-range correlations, nucleon form factors, and higherorder nucleonic weak currents are included in our results. In this work we have also extended the old MCM formalism by adding a previously omitted term in the two-phonon transition density. We also checked the nuclear spectroscopy of the involved daughter nuclei by comparing our nuclear-structure results with the available data on low-energy spectra and electric quadrupole transitions involving the nuclear states of relevance in the present study.
We found in the calculations that the adopted nuclearstructure framework describes reasonably the involved nuclear states and that in the case of the light-neutrino exchange the earlier omission of the added two-phonon term is justified for the heavier nuclei (A = 116, 124, 130, 136) but causes quite significant percentage shifts (up to 44%) for the lighter nuclei. For the heavy-neutrino exchange the added term caused very large shifts in the magnitudes of the NMEs and its omission was no longer realistic.
The decay of 96 Zr to the first 0 + excited state in 96 Mo was found to be amplified relative to the other excited-state decays. This amplification might make this decay transition a promising candidate for future 0νββ experiments.
We have summarized our results in the expressions (27) and (28), connecting the 0νββ half-life measurements with the effective mass of the light Majorana neutrino and the effective inverse mass of the heavy Majorana neutrino via tabulated nuclear-structure coefficients.

ACKNOWLEDGMENTS
This work has been supported partly by the Academy of Finland under the Finnish Centre of Excellence Programme 2012-2017 (Nuclear and Accelerator Based Programme at JYFL) and partly by the Ellen and Artturi Nyyssönen Foundation.