Two-neutrino $\beta\beta$ decay of $^{136}$Xe to the first excited $0^+$ state in $^{136}$Ba

We calculate the nuclear matrix element for the two-neutrino $\beta\beta$ decay of $^{136}$Xe into the first excited $0^+$ state of $^{136}$Ba. We use different many-body methods: the quasiparticle random-phase approximation (QRPA) framework, the nuclear shell model, the interacting boson model (IBM-2), and an effective field theory (EFT) for $\beta$ and $\beta\beta$ decays. While the QRPA suggests a decay rate at the edge of current experimental limits, the shell model points to a half-life about two orders of magnitude longer. The predictions of the IBM-2 and the EFT lie in between, and the latter provides systematic uncertainties at leading order. An analysis of the running sum of the nuclear matrix element indicates that subtle cancellations between the contributions of intermediate states can explain the different theoretical predictions. For the EFT, we also present results for two-neutrino $\beta\beta$ decays to the first excited $0^+$ state in other nuclei.


Introduction
Two-neutrino double-beta (2νββ) decay and twoneutrino double-electron capture (2νECEC) change the atomic number of a nucleus by two via emission or capture of two electrons accompanied by emission of two neutrinos. These second-order weak decays are the rarest processes observed to date, with measured half-lives exceeding 10 21 years [1]. Hence, detecting them demands monitoring tons of otherwise-stable atomic nuclei over several years.
2νββ-decay and 2νECEC half-lives depend on a nuclear matrix element, M 2ν , which encodes information on the structure of the involved nuclei [2]. Therefore, 2νββ-decay measurements can test theoretical predictions of different nuclear many-body calculations of M 2ν . This is valuable because the same methods can also predict the nuclear matrix elements of the neutrinoless double-beta (0νββ) decay, which has not been observed yet. Detecting 0νββ decay promises to unveil the nature of neutrinos, will establish the violation of lepton-number conservation in the laboratory and in general indicate new physics beyond the Standard Model (BSM) of particle physics [3]. Reliable 0νββ nuclear matrix elements are thus needed to extract BSM physics from measurements. However, various many-body methods currently predict matrix elements that differ by up to a factor of a few [3] and almost all calculations miss an important recently acknowledged two-nucleon operator [4][5][6]. Since both 2νββ and 0νββ decays share initial and final states, are sensitive to spin and isospin operators, and may have correlated matrix elements [7,8], tests of M 2ν could help to reduce the theoretical uncertainties of 0νββ nuclear matrix elements.
However, calculating M 2ν is challenging. In fact, seldom has nuclear theory been able to predict 2νββ half-lives before their measurement; calculations yield shorter half-lives than experiment and are corrected with an ad hoc reduction of M 2ν , usually known as "quenching" [9][10][11]. A similar correction is needed in Gamow-Teller (GT) β decays [12][13][14][15], except in ab initio calculations which include two-body currents and many-body correlations [16]. Unfortunately, the latter are not yet capable of reproducing 2νββ decays, even for the lightest emitter 48 Ca [17,18]. For other methods, if the deficiency is systematic [12][13][14], the quenching needed for an unmeasured 2νββ decay can be inferred from other GT and 2νββ decays in the same mass region [10]. Another approach is to construct an effective field theory (EFT) for β decays with low-energy couplings (LECs) adjusted to GT transitions, so that M 2ν is predicted with theoretical uncertainties [19]. These strategies led to the predictions of the 2νββ half-life of 48 Ca by the nuclear shell model [20], and of the 2νECEC rate for 124 Xe by the shell model [21], quasiparticle random-phase approximation (QRPA) [15,22] and EFT [21], in good agreement with subsequent measurements [23,24].
2νββ decays to excited states have already been measured in 100 Mo [1,25] and 150 Nd [1,26]-both without a prior theoretical prediction-and are currently being explored in 76 Ge [27,28], 82 Se [29], 130 Te [30], and 136 Xe [31,32]. In this work, we predict the half-life of the 2νββ decay of 136 Xe to the first 0 + 2 excited state in 136 Ba by several different many-body methods, following similar strategies as for 2νββ decays to the 0 + gs ground state. We provide the first shell-model predictions for this decay using the same Hamiltonians as in previous 136 Xe studies [10,[33][34][35][36]. We also present EFT results with systematic theoretical uncertainties for decays to the 0 + 2 state following Refs. [21,37], thereby also extending the EFT calculations to 136 Xe not included earlier. In addition, we improve previous QRPA results by using larger bases, consistent with 0νββdecay work [36], and we update the interactingboson-model (IBM-2) prediction by using refined model parameters [38]. Figure 1 summarizes our predictions, compared to previous works [11,15] and current experimental limits [31,32]. While our results are consistent with experiment, only a small part of the QRPA band lies inside the non-excluded region. In contrast, the EFT and IBM-2 favor half-lives almost an order of magnitude longer than the current experimental lower limits. Finally, the shell model  [11,15] (in gray) and experimental limits [31,32] (horizontal lines with arrows).
suggests that detecting the decay requires improving the current sensitivities by over two orders of magnitude. These diverse predictions indicate that the 136 Xe 2νββ decay to the 0 + 2 state in 136 Ba will be a very useful test of many-body methods used to calculate 0νββ-decay nuclear matrix elements.

Nuclear matrix element
The 2νββ-decay and 2νECEC nuclear matrix element is given by [2,15] where indices a, b run over all nucleons, the isospin operator τ − turns neutrons into protons, σ is the spin operator, and the denominator involves the electron mass m e and the energies E of the initial (i), final (f ) and each kth intermediate 1 + also depends on a well-known phase-space factor G 2ν [41,42] and the axial nucleon coupling g A . For 136 Xe, the phase space disfavors decays to the 0 + From a measured half-life one can obtain the socalled effective matrix element [1]: which carries all the information to be extracted from nuclear theory. However, our shell-model, QRPA and IBM-2 results require quenching. In these cases, we include a quenching factor q obtained matching the calculated M 2ν to the experimentally recommended M 2ν eff value [1] for groundstate-to-ground-state decay, so that we effectively

Quasiparticle random-phase approximation
The QRPA method considers the 0 + ground states of the initial and final ββ nuclei as QRPA vacua, |QRPA , building nuclear excitations on top of them. The ββ-decay nuclear matrix elements are computed summing over intermediate states in the odd-odd nucleus [see Eq. (1)], which are obtained by performing proton-neutron (pn) QRPA diagonalizations based on the initial and final states as where J π k denotes the spin-parity of the kth intermediate state with omitted projection quantum number M , a † ( a) are nucleon creation (annihilation) operators, and X and Y the pnQRPA backward and forward amplitudes.
On the other hand, the first excited 0 + 2 state in 136 Ba can be described in the charge-conserving (cc) QRPA formalism as a two-phonon excitation: where with normalization factor N ab = 1/ √ 1 + δ ab . We follow Ref. [43] to obtain the transition densities needed to evaluate the matrix element in Eq. (1) using the multiple-commutator model. This formalism includes a correction to the two-phonon transition density that was omitted in an earlier QRPA study [15].
The pnQRPA is adjusted in the usual way: the particle-hole parameter g ph reproduces the energy of the GT giant resonance in 136 Cs [44], and the particle-particle parameter g pp is fitted to the 2νββ-decay half-life to the 0 + gs of 136 Ba, following the partial isospin-restoration scheme [45]. We assume a quenching range q = (0.47 − 1.00), covering typical values of the effective axial coupling g eff A = qg A = (0.6 − 1.27) [15]. For the decay to an excited state there are two additional ccQRPA parameters. Due to its nature, the 2 + 1 state is mainly sensitive to the particle-hole parameter G ph , which we adjust so that the calculated 2 + 1 energy equals half of the measured energy of the 0 + 2 , 2 + 2 , 4 + 1 twophonon triplet. Since the three energies are not exactly degenerate, we obtain a range of values for G ph . For the particle-particle parameter we keep the default value G pp = 1.0 [15,43] because it does not affect much the 2 + 1 state. We consider large no-core single-particle bases consisting of 26 single-particle orbitals [36,44], using a Woods-Saxon and also an adjusted basis. The latter has an increased spin-orbit splitting of the neutron 0h orbitals by 1.5 MeV to better reproduce excitation energies in the neighboring odd-mass nuclei.
The low-lying energy spectra of 136 Xe obtained with the GCN5082 interaction agrees very well with experiment [58], and for QX energies are a little too high [by ∼ 200 keV for the 2 + 1 state, ∼ 300 keV for the 4 + 1 and 6 + 1 states] but also of good quality. 3 GCN5082 also gives a better overall description for 136 Ba, but the excitation energy of the 0 + 2 state agrees with data in a similar way for both interactions: 1.44 MeV (GCN5082) and 1.80 MeV (QX) compared to the experimental 1.58 MeV.
For the M 2ν calculations we use Eq. (1) through the explicit computation of 1 + states in 136 Cs using the strength function method [46]. As it is well known, the nuclear shell model overpredicts M 2ν for the decay to the 0 + gs , and quenching is needed to reproduce experiment: q = 0.42 (GCN2850) [10] and q = 0.68 [33,35] (QX). The same quenching factors are used for the decay to the 0 + 2 state.

Microscopic interacting boson model
The IBM-2 [59,60] maps the fermion Hamiltonian onto a boson space [61] and evaluates it with realistic bosonic wave functions. The method is discussed in detail in previous ββ-decay studies [11,62], which include a calculation of the 136 Xe 2νββ decay to the 0 + 2 state in 136 Ba [11]. Here we improve the IBM-2 calculation using reassessed single-particle and -hole energies and interaction strengths which lead to single-particle occupancies in better agreement with nucleon-removal experiments [38]. We fit the IBM-2 parameters to reproduce the spectroscopic data of the low-lying energy states for 136 Xe, and for 136 Ba we take the parameters from Ref. [63]. The IBM-2 configuration space is the same as for the shell model: 1d 5/2 , 0g 7/2 , 2s 1/2 , 1d 3/2 , and 0h 11/2 for both neutrons and protons. One should note that the IBM-2 matrix elements to the 0 + 2 state are sensitive to the so-called IBM-2 Majorana parameters: In Ref. [64] it was found that the corresponding 150 Nd 0νββdecay matrix element doubled when adjusting them to new data on 150 Sm scissors-mode transitions to excited 0 + and 2 + states. Unfortunately, similar measurements for 136 Ba are not available.
The IBM-2 calculations assume the closure approximation, replacing the energies of the intermediate states in Eq. (1) with an average energy E k , and then summing analytically over intermediate states. Thus, Eq. (1) is simplified to where E k = 11.32 MeV for 136 Cs. Like in the nuclear shell model, the 2νββ-decay matrix element obtained using Eq. (7) is overestimated, and a quenching factor q = 0.31 is used to agree with the  Figure 2: EFT 2νββ nuclear matrix elements (2νECEC for 124 Xe) with theoretical uncertainties at leading order, obtained with LECs fitted to GT data as in Ref. [21] (hatched blue bars) [67][68][69][70][71][72][73][74] and to 2νββ half-lives to the 0 + gs , the strategy of this work (solid red bars) [1,24]. The nuclear matrix elements are compared to empirical values (black circles) [1,24]. Bottom: 0 + gs to 0 + gs decays. Top: 0 + gs to 0 + 2 decays.
empirical value [1]. We assume the same q for the decay to the 0 + 2 state. In addition, for this decay we estimate the sensitivity to parameter changes and model assumptions, including quenching and the closure approximation, to be ±21%, as discussed in detail in Ref. [65].

Effective field theory
The EFT is formulated in terms of nucleon and phonon degrees of freedom coupled to a spherical even-even core. Once the LECs are fitted to experiment (e.g., to β decays), the EFT predicts processes dictated by the same operator (e.g., 2νββ decay). In addition to nuclear properties, the EFT systematically provides the associated theoretical uncertainties based on its power counting [66]. The EFT describes well GT transitions to excited states and 2νββ decays [19], including the prediction of the 2νECEC half-life of 124 Xe [21]. Recently it has been applied to 0νββ decay as well [37].
In the EFT, M 2ν can be calculated using the single-state-dominance approximation [19]. Previous EFT studies fitted the LECs to β decay or GT strengths [19,21], which for 136 Xe are limited to the GT strength from the 136 Xe( 3 He,t) 136 Cs chargeexchange reaction [73]. We label this matching EFT GT . The bottom panel of Fig. 2 (blue bands) shows that, even though this strategy works very well for most of the measured decays within uncertainties [19,21], it overpredicts the 136 Xe 2νββ nuclear matrix element to the 0 + gs . We attribute this deficiency to the fact that the low-energy spectrum of 136 Xe does not exhibit the properties of a spherical collective system [39]. In contrast, 136 Ba shows a one-phonon excited state at approximately half of the energy of the two-phonon excitations (∼ 818/1565 keV) and typical ps lifetimes for the collective excitations [39]. However this nucleus cannot be used to fit the LECs due to lack of GT data. For 48 Ca and 96 Zr the EFT also better describes the final nuclei, but in these cases there is GT information available involving 48 Ti and 96 Mo to fit the LECs. Figure 2 shows that for these decays the EFT agrees well with experiment, but the agreement is lost (for 48 Ca) or worsens (for 96 Zr), if the LECs are only fitted to GT data from the initial nuclei, as in 136 Xe. Finally, the EFT does not describe well the 150 Nd 2νββ decay to the 0 + gs , probably due to its deformation.
Therefore, we here follow an alternative strategy, labelled EFT 2ν , in contrast to the previous EFT GT . In the EFT with single-state-dominance approximation the matrix element for the decays into the 0 + 2 state and 0 + gs can be related by [19] M 2ν The analytical expression of the associated uncertainty is given by Eq. (44) in Ref. [19]. We then fit the LECs directly to measured 2νββ decays effectively through its M 2ν EFT matrix element, assigning an EFT uncertainty that includes the single-state-dominance approximation and the uncertainties in M 2ν EFT from the EFT truncation and the experimental values. By construction, the EFT 2ν reproduces the empirical data for the decay to the 0 + gs , see bottom panel in Fig. 2 (red bands).
Then, we use Eq. (8) to predict the decays to the 0 + 2 excited state in the EFT. The top panel in Fig. 2 shows that this approach agrees well with the two measured decays in 100 Mo and 150 Nd. Therefore, we use the EFT 2ν with its leading order uncertainties as recommended values in this work for 136 Xe. Note that for the decay to the 0 + 2 state the EFT 2ν . We give ranges for the QRPA with Woods-Saxon (WS) and adjusted (adj.) bases, the nuclear shell model (NSM) with the GCN5082 and QX interactions, our IBM-2 results, the EFT with its leading-order uncertainties, as well as previous QRPA and IBM-2 results.

Method
Reference T  Table 1 gives the predicted half-life of the 136 Xe 2νββ-decay to the 0 + 2 state of 136 Ba obtained with the QRPA, nuclear shell model, IBM-2, and EFT. Figure 1 shows the combined ranges for each manybody method. All of them have been adjusted, either via LECs or a quenching factor, to reproduce the known 2νββ decay to the 0 + gs . For the QRPA, the uncertainty range for each basis is dominated by the G ph and quenching ranges considered. In the shell model, the band for each Hamiltonian is simply given by the empirical M 2ν eff [1] used to obtain q. The estimated IBM-2 and calculated EFT 2ν uncertainties are explained in Secs. 3.3 and 3.4, respectively.

Results and discussion
The dispersion of the theoretical predictions is striking: except for the IBM-2 and EFT 2ν -which are fully consistent within uncertainties-different many-body methods predict inconsistent half lives ranging more than four orders of magnitude. Thus, in the extremes, the shell-model and QRPA matrix elements can disagree by up to a factor 100. This difference is much more pronounced than the one for 0νββ-decay matrix elements [3], which does not exceed a factor five or so.
In order to understand these results, Figs and shell model, respectively. References [35,75] performed a similar analysis limited to the decay to the 0 + gs . Figure 3 shows that for the QRPA, for the decay to the 0 + 2 state all intermediate states contribute with the same sign, contrary to the decay to the 0 + gs where there are some cancellations at high energies E ∼ 10 MeV, as previously pointed out in Ref. [75]. In addition, the individual contributions to the decay to the excited state are also larger, especially in the case of the adjusted basis; here the contribution of 1 + states around 2 MeV exceeds the summed contribution of all intermediate states for the ground-state decay. Therefore, the QRPA matrix element is much larger for the decay to the 0 + 2 state, leading to a relatively short half-life. This is especially the case for the adjusted basis, whose prediction is in strong tension with experimental limits, see Table 1 and Fig. 1.
In contrast, Fig. 4 shows that for the shell model there are no cancellations between contributions to M 2ν eff for the decay to the 0 + gs , while intermediate states add up with different signs in the decay to the 0 + 2 state. The qualitative behaviour of the two running sums is similar for the two Hamiltonians. As a consequence, the shell-model matrix element is significantly smaller for the decay to the 0 + 2 state, and the corresponding half-life becomes very long. Note that while Fig. 3   states with energies E ∼ 15 MeV, the shell-model matrix elements in Fig. 4 are converged at E ∼ 10 MeV. This suggests that the shell model could be missing contributions from high-energy states due to its limited configuration space. If these potential contributions had positive sign, like in the QRPA, they would shorten the half-life of the decay to the 0 + 2 state predicted by the shell model. On the other hand, the EFT with single-state dominance does not allow for strong cancellations in the matrix elements, since possible cancellations are just part of the theoretical uncertainty. This is in contrast with our analysis of the QRPA (decay to the 0 + gs ) or shell-model (decay to the 0 + 2 state) running sums in Figs. 3 and 4. This qualitative difference in the running sums cannot be assessed in the lower-resolution EFT using the single-state dominance approximation. For instance, high-energy 1 + states which are relevant for the cancellations of the QRPA and NSM running sums may not be fully captured by the EFT. Nonetheless, the agreement with measured decays to 0 + 2 states shown in Fig. 2 suggests that the EFT theoretical uncertainties at least in part capture the uncertainty associated with this approximation. Due to absence of explicit cancellations, the EFT decay to the 0 + gs and 0 + 2 states can be expected to be relatively similar. The IBM-2 matrix elements for the two decay branches obtained in the closure approximation are also similar. This naturally leads to the EFT and IBM-2 intermediate half-lives in Table 1 and Fig. 1,   6 about 10 3 − 10 4 times longer than the decay to the 0 + gs as dictated by the different phase-space factors.

Summary
In summary, we have predicted the half-life of the 136 Xe 2νββ decay to the 0 + 2 state of 136 Ba using four different many-body methods that are also used to calculate 0νββ-decay nuclear matrix elements. Our results indicate a large uncertainty from nuclear theory: while the QRPA prediction is close to current limits, the nuclear shell model indicates a half-life more than two orders of magnitude longer. The IBM-2 and EFT results lie in between. We have provided error estimates for all four methods, but only the EFT ones can be considered as systematic theoretical uncertainties. Our findings thus highlight that further experimental searches of this decay are very useful tests of theoretical models aiming to predict 0νββ nuclear matrix elements. Finally, for the EFT we have also presented results for 2νββ decays to the excited 0 + 2 state in other nuclei, which agree well within uncertainties for the two measured cases.