Precision Analysis of the 136Xe Two-Neutrino ββ Spectrum in KamLAND-Zen and Its Impact on the Quenching of Nuclear Matrix Elements

KamLAND-Zen Collaboration; Gando, A.; Gando, Y.; Hachiya, T.; Minh, M. Ha; Hayashida, S.; Honda, Y.; Hosokawa, K.; Ikeda, H.; Inoue, K.; Ishidoshiro, K.; Kamei, Y.; Kamizawa, K.; Kinoshita, T.; Koga, M.; Matsuda, S.; Mitsui, T.; Nakamura, K.; Ono, A.; Ota, N.; Otsuka, S.; Ozaki, H.; Shibukawa, Y.; Shimizu, I.; Shirahata, Y.; Shirai, J.; Sato, T.; Soma, K.; Suzuki, A.; Takeuchi, A.; Tamae, K.; Ueshima, K.; Watanabe, H.; Chernyak, D.; Kozlov, A.; Obara, S.; Yoshida, S.; Takemoto, Y.; Umehara, S.; Fushimi, K.; Hirata, S.; Berger, B.E.; Fujikawa, B.K.; Learned, J. G.; Maricic, J.; Winslow, L.A.; Efremenko, Y.; Karwowski, H.J.; Markoff, D. M.; Tornow, W.; O'Donnell, T.; Detwiler, J.A.; Enomoto, S.; Decowski, M.P.; Menéndez, J.; Dvornický, R.; Šimkovic, F.

We present a precision analysis of the 136 Xe two-neutrino ββ electron spectrum above 0.8 MeV, based on highstatistics data obtained with the KamLAND-Zen experiment. An improved formalism for the two-neutrino ββ rate allows us to measure the ratio of the leading and subleading 2νββ nuclear matrix elements (NMEs), ξ 2ν 31 ¼ −0.26 þ0.31 −0. 25 . Theoretical predictions from the nuclear shell model and the majority of the quasiparticle random-phase approximation (QRPA) calculations are consistent with the experimental limit. However, part of the ξ 2ν 31 range allowed by the QRPA is excluded by the present measurement at the 90% confidence level. Our analysis reveals that predicted ξ 2ν 31 values are sensitive to the quenching of NMEs and the competing contributions from low-and high-energy states in the intermediate nucleus. Because these aspects are also at play in neutrinoless ββ decay, ξ 2ν 31 provides new insights toward reliable neutrinoless ββ NMEs.
Introduction.-Double-beta (ββ) decay is a rare nuclear process. The ββ decay emitting two electron antineutrinos and two electrons (2νββ) is described within the standard model of the electroweak interaction. In contrast, the ββ mode without neutrino emission (0νββ) implies new physics and can only occur if neutrinos are Majorana particles. While 2νββ decay has been measured in 12 isotopes [1], an observation of 0νββ decay remains elusive. In the standard scenario, the 0νββ rate is proportional to the square of the effective Majorana neutrino mass, m ββ [2], allowing the establishment of definite benchmarks toward the discovery of 0νββ decay in experiments.
The 0νββ rate, however, also depends on nuclear matrix elements (NMEs) which are poorly known [3], as 0νββ NME estimates vary between the many-body approaches used to calculate them. In addition, NMEs may be affected by a possible "quenching" or, equivalently, an effective value of the axial-vector coupling g eff A in the decay. Overall, the NME uncertainty can reduce the experimental sensitivity on m ββ by up to a factor of 5 [4]. To mitigate this, nuclear many-body predictions need to be tested in other observables. Several nuclear structure [5][6][7][8] and Gamow-Teller (GT) properties [9][10][11] have been proposed as 0νββ decay probes. Because 2νββ and 0νββ decays share initial and final nuclear states, and the transition operators are similar, a reproduction of 2νββ decay is key to reliable 0νββ NME predictions. Nonetheless, few nuclear manybody methods are well suited for both ββ modes, because nuclei with even and odd numbers of neutrons and protons up to high excitation energies need to be described consistently. The most notable approaches are the quasiparticle random-phase approximation (QRPA) [12][13][14][15][16] and the nuclear shell model [17][18][19][20][21][22].
The 2νββ rate is usually expressed as where M 2ν GT is the 2νββ NME and G 2ν 0 a known phase-space factor [23]. As a result, g eff A can be determined from the measured T 2ν 1=2 once M 2ν GT is theoretically evaluated, a strategy followed in Ref. [24]. While a similar approach has been used in the nuclear shell model, especially for 136 Xe [20,25], it is more common to take g eff A from GT β decay and electron-capture (EC) rates [25,26], assuming a common quenching for all weak processes. Likewise, the QRPA can also use β decay and EC to obtain g eff A [27][28][29], even though the standard approach is to fix g eff A first, and then adjust the nuclear interaction so that M 2ν GT describes the 2νββ half-life [30]. In this way, the nuclear shell model and QRPA typically reproduce experimental 2νββ rates and predict nonmeasured ones [17,18,[31][32][33].
Recently, the 2νββ decay of several isotopes has been observed with high statistics by the NEMO-3 [34], EXO [35], KamLAND-Zen [36], GERDA [37], MAJORANA [38], and CUORE [39] collaborations. These achievements demand an improved theoretical description. Reference [40] gives a more accurate expression for the 2νββ decay rate where the phase-space factor G 2ν 2 has a different dependence on lepton energies than G 2ν 0 , and the subleading nuclear matrix element M 2ν GT−3 enters the (real-valued) ratio ξ 2ν GT is sensitive to contributions from high-lying states in the intermediate odd-odd nucleus, for M GT−3 only the lowest-energy states are relevant due to rapid suppression in the energy denominator. Consequently, ξ 2ν 31 probes additional, complementary physics to the 2νββ half-life. This novel observable can be determined experimentally by fitting the 2νββ electron energy spectrum to extract the leading and second order contributions in Eq. (2). Hence, the measurement of ξ 2ν 31 challenges theoretical calculations and can discriminate between those that reproduce the 2νββ rate.
In this Letter, we analyze the high-statistics 2νββ decay of 136 Xe with KamLAND-Zen [36] and compare the measured T 2ν 1=2 and ξ 2ν 31 values with the predictions from the QRPA and nuclear shell model. In KamLAND-Zen, the spectral distortion due to ξ 2ν 31 could be up to 8% based on the theoretical predictions. Such effect is testable with accumulated statistics of ∼10 5 2νββ decays. Because 0νββ NMEs also show a competition between contributions from low-and high-energy intermediate states [15], testing theoretical ξ 2ν 31 predictions can provide new insights on 0νββ calculations, including the possible quenching of the NMEs.
Experiment and results.-The KamLAND-Zen (KamLAND Zero-Neutrino Double-Beta Decay) detector consists of 13 tons of Xe-loaded liquid scintillator (Xe-LS) contained in a 3.08-m-diameter spherical inner balloon (IB). The IB is constructed from 25-μm-thick transparent nylon film and is suspended at the center of the KamLAND detector [41,42]. The IB is surrounded by 1 kton of liquid scintillator (LS) which acts as an active shield. The scintillation photons are viewed by 1879 photomultiplier tubes mounted on the inner surface of the containment vessel. The Xe-LS consists of 80.7% decane and 19.3% pseudocumene (1,2,4-trimethylbenzene) by volume, 2.29 g=L of the fluor PPO (2,5-diphenyloxazole), and ð2.91 AE 0.04Þ% by weight of enriched xenon gas. The isotopic abundances in the enriched xenon were measured by a residual gas analyzer to be ð90.77 AE 0.08Þ% 136 Xe, ð8.96 AE 0.02Þ% 134 Xe.
We report on data collected between December 11, 2013 and October 27, 2015, which is the same data set analyzed for the 0νββ search in Ref. [36] with a total live time of 534.5 days. The selection to reduce the background contributions is the same as in Ref. [36], but we apply a tightened 2νββ event selection for this work in order to avoid systematic uncertainties arising from backgrounds. The fiducial volume for the reconstructed event vertices is defined as a 1-m-radius spherical shape at the detector center, which gives a fiducial exposure for this analysis of ð126.3 AE 3.9Þ kg yr in 136 Xe. We perform a likelihood fit to the binned energy spectrum of the selected candidates between 0.8 and 4.8 MeV, tightened relative to the 2νββ analysis in Ref. [36]. The systematic uncertainties on the 2νββ rate are evaluated identically as in Ref. [36] and are summarized in Table I. A detailed energy calibration is essential for the extraction of ξ 2ν 31 . The energy scale was determined using γ rays from 60 Co, 68 Ge, and 137 Cs radioactive sources, γ rays from the capture of spallation neutrons on protons and 12 C, and β þ γ-ray emissions from 214 Bi, a daughter of 222 Rn (lifetime 5.5 day) that was introduced during the Xe-LS purification. Uncertainties from the nonlinear energy response due to scintillator quenching and Cherenkov light production are constrained by the calibrations. The most important calibration is the high-statistics 214 Bi from the initial 222 Rn distributed uniformly over the Xe-LS volume. To ensure that the calibration with 214 Bi can be applied to the entire data set, we confirmed that the time variation of the energy scale is less than 0.5% based on the spectral fit to the 2νββ decays for each time period. This uncertainty is reduced relative to the previous analysis [36], and is added to the energy scale error, which is the dominant error source for the ξ 2ν 31 measurement, as discussed later. The energy spectrum of selected candidate events between 0.8 and 2.5 MeV together with the best-fit spectral decomposition is shown in Fig. 1. In the fit, the contributions from 2νββ and major backgrounds in the Xe-LS, such as 40 K, 210 Bi, and the 228 Th-208 Pb subchain of the 232 Th series are free parameters and are left unconstrained. The background contribution from 110m Ag, which is important for the 0νββ analysis, is also a free parameter in the fit. The contributions from the 222 Rn-210 Pb subchain of the 238 U series, and from 11 C and 10 C (muon spallation products), as well as the detector energy response model parameters, are allowed to vary but are constrained by their independent estimations [36].
The 2νββ spectrum is computed with Eq. (2), convolved with the detector response function. It is characterized by two free parameters: the total 2νββ rate and the ratio of the matrix elements ξ 2ν 31 . We obtained a best fit of ξ 2ν 31 ¼ −0.26 þ0. 31 −0.25 and a 90% confidence level (C.L.) upper limit of ξ 2ν 31 < 0.26. The systematic uncertainty on the energy scale limits the sensitivity of the ξ 2ν 31 measurement, because an energy scale shift introduces a shape distortion similar to the change generated by a nonzero ξ 2ν 31 . The best-fit total 2νββ rate in the Xe-LS mass is 99.7 þ1.2 −1.4 ðton dayÞ −1 . Figure 2 shows the joint confidence intervals for the 2νββ rate and ξ 2ν 31 , which exhibit only a slight positive correlation. It indicates that the effect on the total 2νββ rate estimate by the introduction of the second order contribution is small. The effect on the 0νββ analysis is also negligibly small. Considering the systematic uncertainties in Table I   T 2ν 1=2 ¼ 2.165 AE 0.016ðstatÞ AE 0.059ðsystÞ × 10 21 yr [35]. Our analysis neglects counts from decays to excited states in 136 Ba, for which our shell model calculations predict T 2ν 1=2 ð0 þ gs → 0 þ 1 Þ > 10 26 yr. Even a very conservative halflife of 8.7 × 10 24 yr that assumes the same NME for the decay to the ground (gs) and excited 0 þ states does not affect our results. This issue might need to be revisited in the case of an unexpectedly short half-life close to the present 90% C.L. lower limit of 8.3 × 10 23 yr [43]. The correction to 2νββ decay represented by ξ 2ν 31 impacts KamLAND-Zen analyses of spectral distortions, including extraction of half-lives to excited states as well as searches for beyond-standard-model physics, such as for Majoron emission modes. Considering ξ 2ν 31 as a free parameter, we find the additional uncertainty comparable to the energy scale error. Updated spectral analyses will be presented in future publications.
with energy denominator Δ ¼ ½E j − ðE i þ E f Þ=2=m e . E k is the energy of the nuclear state jJ π k i with total angular momentum J and parity π, and m e is the electron mass. The labels i, j, f refer to the initial, intermediate, and final nuclear states, respectively, while σ is the spin and τ − the isospin lowering operator.
We perform nuclear shell model calculations in the configuration space comprising the 0g 7=2 , 1d 5=2 , 1d 3=2 , 2s 1=2 , and 0h 11=2 single-particle orbitals for both neutrons and protons, using the shell model code NATHAN [44]. We reproduce M 2ν GT ¼ 0.064 from Ref. [25] with the GCN interaction [19] and also use the alternative MC interaction from Ref. [45], which yields M 2ν GT ¼ 0.024. Both interactions have been used in 0νββ decay studies [11,46]. Shell model NMEs for β and 2νββ decays are typically too large, due to a combination of missing correlations beyond the configuration space, and neglected two-body currents in the transition operator [3]. This is phenomenologically corrected with a "quenching" factor q, or g eff A ¼ qg A . In general, the quenching that fits GT β decays and ECs in the same mass region is valid for 2νββ decays as well. Around 136 Xe, GT transitions with GCN are best fit with q ¼ 0.57 [25], and with the same adjustment the 136 Xe GT strength into 136 Cs [10], available up to energy E ≲ 4.5 MeV, is well reproduced by both interactions. However, the experimental 2νββ halflife suggests different quenching factors q ¼ 0.42ð0.68Þ for GCN (MC). The calculations yield M 2ν GT−3 ¼ 0.011ð0.0025Þ. We assume a common quenching for M 2ν GT and M 2ν GT−3 because the shell model reproduces well GT strengths at low and high energies up to the GT resonance [9]. This gives ratios ξ 2ν 31 ¼ 0.17 for GCN and ξ 2ν 31 ¼ 0.10 for MC, both consistent with the present experimental analysis.
We also perform 2νββ decay QRPA calculations with partial restoration of isospin symmetry [16]. We consider a configuration space of 23 single-particle orbitals (the six lowest harmonic oscillator shells with the addition of the 0i 13=2 and 0i 11=2 orbitals). We take as nuclear interactions two different G matrices, based on the charge-dependent Bonn (CD-Bonn) and the Argonne V18 nucleon-nucleon potentials. We fix the isovector proton-neutron interaction imposing the restoration of isospin [16]. Finally, we adjust the isoscalar neutron-proton interaction to reproduce the 2νββ decay half-life for different values in the range Discussion.- Figure 3 shows the effective axial-vector coupling constant g eff A as a function of the matrix element

M 2ν
GT−3 for the 2νββ decay of 136 Xe. A large region in the g eff A − M 2ν GT−3 plane is excluded by the present 90% C.L. limit ξ 2ν 31 < 0.26. The two nuclear shell model GCN and MC results, indicated by points, are consistent with the KamLAND-Zen limit. The QRPA Argonne and CD-Bonn results are presented by curves, which accommodate 0.33 ≤ g eff A ≤ 1.269 values (the lower end corresponds to vanishing isoscalar interactions). Both curves are very similar, because QRPA ratios of matrix elements with the same initial and final states are weakly sensitive to the nucleon-nucleon interaction [30]. Figure 3 shows that, even though most QRPA predictions are consistent with our measurement, g eff A ≳ 1.14ð1.00Þ for the Argonne (CD-Bonn) potential is excluded at 90% C.L. by the KamLAND-Zen ξ 2ν 31 limit. Figure 3 also shows that for g eff A ≳ 0.7 the QRPA predicts larger ξ 2ν 31 values than the nuclear shell model. Elsewhere, the QRPA ratios lie between those of the GCN and MC shell model interactions. Interestingly, for g eff A ∼ 0.5, the QRPA and shell model GCN results are close. While such relatively small g eff A values are not always considered in 2νββ QRPA calculations of 136 Xe, they are favored by QRPA statistical analyses that take into account experimental EC and β rates [27,47].
To illustrate the origin of the differences between the theoretical calculations, Fig. 4  GT exceeds the shell model significantly, with a strong g eff A sensitivity. While for g eff A ¼ 1.269 the maximum of the QRPA running sum is almost four times larger than the shell model one, for g eff A ∼ 0.5-not shown in Fig. 4-the difference is only about 20%, consistent with the more similar ξ 2ν 31 values predicted. E exc ∼ 10 MeV shellmodel contributions may be too small due to missing spinorbit partner orbitals, but the QRPA may also overestimate them. Measurements of charge-exchange reactions up to the 136 Xe GT resonance, currently limited to lower energy [10,48], can clarify this picture. Above E exc ≳ 10 MeV, the QRPA excess with respect to the shell model is canceled. The final value, set by the 2νββ half-life, is common to all calculations.
By contrast, Fig. 4 shows that in both shell model and QRPA the lowest 1 þ state component dominates the M 2ν GT−3 NME. Such contribution is more salient for the shell model GCN and QRPA g eff A ¼ 1.269 calculations, which explains the larger associated ξ 2ν 31 value compared to the shell model MC and QRPA g eff A ¼ 0.8 results, respectively. The contrast in the M 2ν GT−3 running sum at low energies is ultimately responsible for the different ξ 2ν 31 values predicted by the QRPA and nuclear shell model.   In 0νββ decay, the running sum of the NME can extend to even higher energies, because in this case there is no dependence on the energy of the intermediate states in the denominator; see Eqs. (3) and (4). Therefore, a competition between contributions from low-and highenergy states similar to 2νββ decay is expected [15,49,50]. Consequently, fixing ξ 2ν 31 in 2νββ decay will allow one to identify the most promising 0νββ NME predictions.
Further experimental ξ 2ν 31 sensitivity improvements may distinguish between various scenarios. On the one hand, measured values of ξ 2ν 31 ≥ 0.11 will allow QRPA calculations to fix the quenched value of g eff A , reducing uncertainties in QRPA 0νββ NMEs [30,51,52]. Likewise, a measured value ξ 2ν 31 ≃ 0.17ð0.10Þ would suggest that the GCN (MC) shell model interaction, with its associated g eff A value, leads to a more reliable 0νββ NME. Because the QRPA and shell model rely on different assumptions, and for 2νββ decay they can exhibit contrasting sensitivities on g eff A -as shown in Fig. 4-a measurement of ξ 2ν 31 could lead to different g eff A values for each model. Furthermore, the quenching may not be the same in 2νββ and 0νββ decays, especially in the light of the differences in the two-body [53][54][55] and contact [56] corrections to the two ββ transition operators. On the other hand, a small ratio ξ 2ν 31 < 0.11, which cannot be accommodated in the present QRPA calculations, or a determination of ξ 2ν 31 very different to the GCN and MC predictions, would demand improved theoretical developments.
Summary.-We have presented a precision analysis of the 136 Xe 2νββ electron spectrum shape with the KamLAND-Zen experiment. For the first time, we set a limit on the ratio of nuclear matrix elements ξ 2ν 31 < 0.26 (90% C.L.). The experimental limit is consistent with the predictions from the nuclear shell model and most QRPA calculations, but excludes QRPA Argonne (CD-Bonn) results for g eff A ≳ 1.14(1.00). The allowed theoretical values vary in the range ξ 2ν 31 ¼ ð0.10 − 0.26Þ, so that future ξ 2ν 31 measurements will be required to further test 2νββ calculations, and select the most successful ones. The associated g eff A value, or NME quenching, would also be identified. Future experiments such as KamLAND2-Zen [57] and others with improved resolution and reduced backgrounds promise enhanced sensitivity to reach this goal. Our analysis reveals that ξ 2ν 31 is sensitive to competing contributions to the NME from low-and high-energy intermediate states. Because a similar competition is also relevant for 0νββ decay, studies of this observable provide new insights for identifying reliable 0νββ NMEs.