Benchmarking $^{136}$Xe Neutrinoless $\beta\beta$ Decay Matrix Element Calculations with the $^{138}{\rm Ba}(p,t)$ Reaction

We used a high-resolution magnetic spectrograph to study neutron pair-correlated $0^+$ states in $^{136}$Ba, produced via the $^{138}{\rm Ba}(p,t)$ reaction. In conjunction with state-of-the-art shell model calculations, these data benchmark part of the dominant Gamow-Teller component of the nuclear matrix element (NME) for $^{136}$Xe neutrinoless double beta ($0\nu\beta\beta$) decay. We demonstrate for the first time an evaluation of part of a $0\nu\beta\beta$ decay NME by use of an experimental observable, presenting a new avenue of approach for more accurate calculations of $0\nu\beta\beta$ decay matrix elements.

The massive nature of neutrinos leads to a violation of the γ 5 invariance [1] for weak interactions. Consequently, there is substantial interest worldwide [2][3][4] to search for standard-model-forbidden neutrinoless double beta (0νββ) decays, that violate lepton number conservation by 2 units. The observation of such decays would prove that the electron neutrino (ν e ) is a Majorana fermion, and therefore indistinguishable from its antiparticle (ν e ). This is consistent with most theories beyond the standard model [5], that attribute the smallness of neutrino masses to a violation of total lepton number at an energy scale of ∼ 10 15 GeV [3,5].
If the mechanism driving a 0νββ decay is via the exchange of a light left-handed Majorana neutrino, then the decay amplitude is proportional to where m ee is the effective Majorana mass of the electron neutrino and M 0ν is the nuclear matrix element (NME) for the decay. The NME is expressed as the sum of Gamow-Teller (GT), Fermi (F) and Tensor (T) components where the Gamow-Teller contribution is the dominant term. In Eq. (1), the U ej are elements of the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) neutrino mixing matrix [6,7], the m j 's are the light neutrino masses and the α j 's are phases in the mixing matrix. For the special case of three-neutrino mixing, the PMNS matrix is parameterized in terms of three mixing angles and one Dirac and two Majorana CP-violating phases [2]. It is evident from Eq. (1) that in addition to the observation of a 0νββ decay process, it is also equally important to determine its half-life, which would establish the absolute neutrino mass scale. Furthermore, such a measurement also has the potential to identify the correct neutrino mass spectrum [8] and find extra sources of CP-violation in the leptonic sector [9]. However, achieving the above (or placing stringent constraints on any new physics) requires an accurate evaluation of the NME for the decay. This has been at the forefront of nuclear physics research in recent times, with several approaches being used to calculate the NMEs for 0νββ decay candidates [10,11]. Depending on the method used, the calculations for specific isotopes disagree with one another, differing by factors of three or more in many cases [10,11]. These discrepancies result  [21] 2.91 Deformed Skyrme-HFB-QRPA (Chapel Hill) [22] 1.55 Spherical QRPA (Bratislava-Tübingen-CalTech) [23] 2.46 ISM (Strasbourg-Madrid) [24] 2.19 ISM (Michigan) [12] 1.46 CDFT [25] 4.24 NREDF [26] 4.77 IBM-2 [27] 3.05 GCM [28] 2.35 in large uncertainties for the NMEs, which not only limit the physics that can be addressed, but also the planning and execution of future 0νββ decay experiments [11]. In contrast, the NMEs for the rare (yet standard-modelallowed) two-neutrino double beta (2νββ) decays can be extracted directly from measured half-lives. These and other experimentally derived spectroscopic information have played a critical role in constraining various NME calculations [12][13][14][15][16][17].
One of the most promising candidates for observing 0νββ decays is 136 Xe ββ → 136 Ba. Its 2νββ decay half-life is much longer than most other cases [18]. As a result, the ratio of the 0νββ decay signal to the irreducible 2νββ decay background in the vicinity of the decay endpoint energy is expected to be larger for this particular case. In fact, a highly sensitive search for 0νββ decays was recently reported for 136 Xe by the KamLAND-Zen collaboration [19], who placed the most stringent upper limits to date on the effective neutrino mass m ee < 61 − 165 meV, depending on the choice of NME used.
We list some recent evaluations of M 0ν for 136 Xe ββ decay in Table I. While some of these results are in reasonable agreement with each other, there still exist large discrepancies in the calculated values, depending on the method used. Needless to say, this is a pressing issue as future 136 Xe 0νββ decay experiments aim to improve their sensitivity by at least one order of magnitude [29]. Additionally, next generation experiments also intend to use the method of barium ion tagging [30] in xenon time projection chambers (TPCs). This technique has the potential to reduce room background contributions to insignificant levels, which undoubtedly would further enhance the sensitivity of 136 Xe 0νββ decay experiments over other experimental searches.
Finally, it is expected that due to these advantages, future 136 Xe ββ decay experiments would also place the tightest constraints on possible CP-violating Majorana phases in the PMNS matrix [31]. The phases in the neutrino mixing matrix are potential leptogenesis parameters that may explain the observed baryon asymmetry in the universe [32][33][34].
Quite naturally, due to the several reasons listed above, the 0νββ decay of 136 Xe presents a compelling case to address the accuracy in its calculated NME . Two methods that have been traditionally used to calculate 0νββ decay NMEs are the interacting shell model (ISM) and the quasiparticle random phase approximation (QRPA). Unlike the latter, the shell model calculations use a limited configuration space that is comprised of relatively fewer single-particle states in the vicinity of the Fermi surface. Despite this restriction, large-scale ISM calculations allow arbitrarily complex correlations between the valence nucleons. On the other hand, the QRPA calculations make use of a much larger model space with comparatively simpler configurations. In general, the ISM calculations are known to yield smaller values of M 0ν , compared to the QRPA [2,11]. This discrepancy has been attributed to different approaches in treating the pairing (or seniority) structure of the nuclear wavefunctions [35,36]. For the most part, previous QRPA calculations assumed spherical ground states for the parent and daughter nuclei, wherein the pairing correlations between like nucleons were taken into account using the Bardeen-Cooper-Schrieffer (BCS) approximation [2,11]. It was only recently that deformed QRPA NME calculations were performed for 0νββ decays [20,22], whose results for 136 Xe are listed in Table I. Compared to the spherical QRPA [21], the deformed calculations yield smaller values for the NME, and are in reasonable agreement with the ISM results. The authors of Refs. [20,22] point out that the suppression of the NME in their calculations is mainly due to differences in the pairing content of the initial and final mean fields. Unlike the spherical QRPA, the deformed calculations accounted for the sharp neutron Fermi surface in 136 Xe due to the neutron number N = 82 shell closure. This curtails the overlap between the BCS wavefunctions and leads to a significant reduction in the calculated NMEs [20,22]. The calculations also suggest that the NME can be even more suppressed if the parent and daughter nuclei have different deformations. Such a scenario will either further reduce [20,22] the QRPA overlap factors mentioned above or result in a similar seniority mismatch between the ISM wavefunctions, due to high-seniority [37] components introduced by the deformation. In comparison, the NME calculations using other many-body approaches such as the nonrelativistic energy density functional (NREDF) theory, covariant density functional theory (CDFT), the interacting boson model (IBM-2) or the generator coordinate method (GCM), predict higher values for the NME (Table I). It has been suggested that these values are most likely overestimated, because of the omission of both collective as well as non-collective correlations, depending on the calculation [11,28,38,39].
In light of the above, precise experimental informa- tion elucidating the properties of 136 Xe and 136 Ba nuclei are crucial to benchmark the NME calculations and further reduce their model dependence. Indeed, differences in the valence nucleon occupancies for these nuclei were recently determined using one nucleon transfer reactions [40,41]. Furthermore, the ground state of 136 Ba is not expected to have a nearly spherical structure as 136 Xe or 138 Ba. The even barium isotopes in the N ≤ 82 region are known to be transitional, displaying a structural evolution from spherical to γ-soft behavior with decreasing neutron number [42][43][44]. In this Letter we discuss neutron pairing correlations in 136 Ba, stud- Ref. [47] Ref. [48] Ref. [49] Ref. [50] Ref. [51] FIG. 2. Measured 138 Ba(p, p) angular distribution from this work (expressed in terms of ratio to the Rutherford cross sections) compared to various DWBA predictions based on different global OMP parameters.
ied with the 138 Ba(p, t) reaction. The experiment was performed at the Maier-Leibnitz-Laboratorium (MLL) in Garching, Germany, where a 1.5 µA, 23 MeV proton beam from the MLL tandem accelerator was incident onto a 40µg/cm 2 thick, 99.9% isotopically enriched 138 BaO target, that was evaporated on a 30 µg/cm 2 carbon backing. The light reaction ejectiles were momentum analyzed using the high-resolution Q3D magnetic spectrograph [45,46], whose solid angle acceptance ranged from 2.3 − 14.6 msr during various stages of the experiment. The detector placed at the focal plane of the Q3D spectrograph consisted of two gas proportional counters and one 7-mm-thick plastic scintillator [46]. A cathode strip foil in the second proportional counter provided position information (with a resolution of ∼ 0.1 mm), while the energy losses in the gas counters and the residual energy deposited in the scintillator allowed for particle identification. The integrated beam intensity was determined using a Faraday cup that was placed at 0 • to the beam, and connected to a Brookhaven Instruments Corporation (BIC) current integrator. For our measurements, we obtained triton angular distributions using four magnetic field settings and at ten spectrograph angle settings, ranging from θ lab = 5 • to 50 • . Fig. 1 shows sample triton spectra obtained from this experiment, where we observe states in 136 Ba up to ∼ 4.6 MeV in excitation energy. The energy resolution of the triton peaks were found to be 10 keV. We also took additional 138 Ba(p, p) elastic scattering data over an angular range of θ lab = 15 • to 115 • , in 5 • steps. These data were used to determine both the effective areal density of the 138 Ba target nuclei, as well as the appropriate global optical model potential (OMP) parameters for the incoming 138 Ba + p reaction channel [47][48][49][50][51]. The latter were used in a zero-range distorted wave Born approximation (DWBA) analysis of our data, for which we used the DWUCK4 code [52] with Woods-Saxon potentials. As shown in Fig. 2, based on a comparison of various DWBA calculations with our elastic scattering measurements, we chose the global proton OMP parameters recommended by Varner et al. [51] for the incoming (proton) channel. For the outgoing 136 Ba + t channel we used the OMP parameters provided by Li et al. [53], as they gave the best agreement with our measured triton angular distribution for the ground state in 136 Ba (c.f. Fig. 3). The transfer form factor was calculated assuming a singlestep pick up of a di-neutron in a singlet relative s-state. For the core-2n coupling we used the global OMP from Ref. [54], whose well depth was adjusted to reproduce the binding energy of each neutron [55].
The above approach was used to perform a comprehensive analysis of the angular distributions for all the triton peaks shown in Fig. 1. We defer a detailed discussion on the analysis to a forthcoming article [56]. In this Letter we only focus on the L = 0 angular distributions, which are critical for studying pair-correlated states in even- even nuclei. Our measured angular distributions identify eight new 0 + states in 136 Ba [57,58]. These results are shown in Fig. 3, together with normalized DWBA cross sections. The data show reasonable agreement with DWUCK4 predictions, except for the well-established 0 + 2 state at 1579 keV [57], where the first minimum occurs at approximately twice the predicted value. This is due to an inherent shortcoming of the DWUCK4 calculations. For example, the calculations ignore multi-step processes and interference from different configurations to the pair transfer amplitude, which can alter the shape of the angular distribution. In Table II we list the measured absolute cross sections for these states at θ lab = 5 • , in addition to the L = 0 transfer strengths to the excited 0 + states relative to the ground state, denoted by ǫ i . The latter were determined for each excited state by the product which was obtained by normalizing the DWBA calculations to the data at forward angles (where the DWBA is best satisfied) for each plot in Fig. 3 and then taking the ratio of the normalization factor for the 0 + i state to the ground state normalization factor. This effectively removes the Q-value dependence and other kinematic effects in determining the relative (p, t) strengths.
Our results in Fig. 3 and Table II show that in addition to a number of hitherto unknown 0 + states in 136 Ba, we observe a large fragmentation of the L = 0, (p, t) strength to these states, with ∼ 30% of the ground state strength concentrated at 2.3 and 2.8 MeV. This manifestly indicates a breakdown of the BCS approximation for neutrons in 136 Ba [13]. A similarly large breakdown was not observed in the 128−134 Ba isotopes [44,59,60]. Nevertheless, such a departure from superfluid behavior is not unexpected in a shape transitional region, particularly around closed shell nuclei [13,61]. Therefore, our results indicate that the ground state wavefunctions of 138,136 Ba could be largely dissimilar due to the 'non-spherical' nature of the latter. Additionally, the data also present evidence of a small pairing gap in 136 Ba, that could occur due to vibrational modes in the pairing field [61,62]. While both possibilities cannot be ruled out, the former will have important ramifications for 136 Xe 0νββ decay NME calculations. Previous work showed that the static quadrupole moment of the first 2 + state in 136 Ba could be as large as −0.19 or +0.25 eb [63][64][65], which does not rule out a significantly deformed 136 Ba ground state.
The results from this experiment also allow tests of the nuclear structure models that are used to calculate the NME for 136 Xe 0νββ decay.
We performed one such test using configuration interaction shell model calculations with the NuShellX code [66]. Unlike the DWUCK4 calculations which merely served to identify the observed 0 + states and determine relative strengths, the shell model calculations were used to assess the absolute (p, t) cross sections, distributed over 0 + states in 136 Ba. The calculations used the five-orbital (0g 7/2 , 1d 5/2 , 1d 3/2 , 2s 1/2 , 0h 11/2 ) valence space for protons and neutrons to determine the wavefunctions for the 0 + ground state of 138 Ba and the lowest fifty 0 + states in 136 Ba. NuShellX was then used to calculate two-neutron transfer amplitudes between these states, that served as input for the Fresco [67] coupled-reaction channels code to generate 138 Ba(p, t) angular distribution predictions.
In the Fresco calculations we used the same OMP parameters as the DWUCK4 calculations (for the proton and triton channels) and took into account the coherent sum of both direct and sequential two-step transfer. The single nucleon transfer amplitudes for the two-step part were obtained assuming one intermediate state in 137 Ba for each of the transferred (n, ℓ, j) values. For the 137 Badeuteron coupling, we used the global OMP parameters recommended by An and Cai [68], based on a comparison with 138 Ba(d, d) angular distribution data (similar to Fig. 2) that we obtained independently. We used three different Hamiltonians for the calculations, which were corrected for core-polarization due to configuration mixing with orbitals outside the model space [69]. The first Hamiltonian is from Ref. [70] and is called sn100pn in the NuShellX interaction library [66]. The second Hamiltonian (that we call sn100t) is very similar to sn100pn, except with minor modifications and was used in Ref. [12] to calculate M 0ν for 136 Xe ββ decay, while the third GCN50:82 [71] Hamiltonian was used by the authors of Ref. [36] to calculate the NME for the decay.
In the left panel of Fig. 4 we compare the running sum of our measured (p, t) cross sections at the most forward angle (where the L = 0 strength is concentrated) to the shell model/Fresco calculations. Similar to our experimental results, the calculations predict the cross section to be dominated by the transition to the ground state in nucleon currents [76], we determine the matrix element to be M 0ν GT (J π = 0 + ) = 5.67. Next we evaluated the NME through the J π = 0 + ground state in 134 Xe, both with and without the core-polarization corrections to the TNA described above. These calculations showed that the five-orbital valence space ISM result for the NME required an enhancement factor of ∼ 1.58, due to the required core-polarization corrections. Such an enhancement significantly increases the magnitude of the NME, so that its revised value is M 0ν GT (J π = 0 + ) = 8.96. To make further comparisons, we also performed a largescale spherical QRPA calculation of the NME (with HOC and the CD-Bonn SRC), using the model parameters of Ref. [23] and 28 orbitals for major oscillator shells with N ≤ 6. This resulted in a value M 0ν GT (J π = 0 + ) = 9.63. It is indeed gratifying that the QRPA result (from using a significantly larger number of orbitals) is closer in value to our revised calculation of the NME, compared to the five-orbital model space result. Clearly, previous configuration interaction (ISM) calculations [36] had underestimated the J = 0 component of the Gamow-Teller NME. We recommend improved calculations of this part of NME as well as the canceling M 0ν GT (J > 0) term, perhaps along the lines of a many-body perturbation theory treatment [77], that takes into account physics contributions from beyond the model space. Since the details of the cancellation between the J = 0 and J > 0 contributions to the NME are model dependent, our result also serves to benchmark future calculations of M 0ν GT (J π = 0 + ) along these lines. We note that in order to make the connection with two-nucleon transfer reaction data, it is important that the results of these calculations be expressed in terms of their J π decomposition.
In summary, this work demonstrates for the first time a direct evaluation of part of a 0νββ decay NME using experimental data. In addition to providing a benchmark for future calculations, it also presents a new avenue of approach for evaluating 0νββ decay NMEs more accurately, motivating similar investigations in other candidates. We also report for the first time a large breakdown of the neutron BCS approximation in an even barium nucleus with N ≤ 82. Our observations motivate a reassessment of the neutron pairing approximation in 136 Ba, used in some NME calculations for 136 Xe 0νββ decay and invite further investigations of the shape of 136 Ba.
We are grateful to Ian Thompson for his assistance with the Fresco calculations. This work was partially