Precision measurement of the 136 Xe two-neutrino ββ spectrum in KamLAND-Zen and its impact on the quenching of nuclear matrix elements

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

We present a precision measurement of the 136 Xe two-neutrino ββ electron spectrum above 0.8 MeV, based on high-statistics 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% C.L. 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. Since these aspects are also at play in neutrinoless ββ decay, ξ 2ν 31 provides new insights towards 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 twelve 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 five [4]. To mitigate this, nuclear many-body predictions arXiv:1901.03871v1 [hep-ex] 12 Jan 2019 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. Since 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 many-body 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 randomphase approximation (QRPA) [12][13][14][15][16] and the nuclear shell model [17][18][19][20][21].
The 2νββ rate is usually expressed as where M 2ν GT is the 2νββ NME and G 2ν 0 a known phasespace factor [22]. 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. [23]. While a similar approach has been used in the nuclear shell model, especially for 136 Xe [20,24], it is more common to take g eff A from GT βdecay and electron-capture (EC) rates [24,25], assuming a common quenching for all weak processes. Likewise, the QRPA can also use β-decay and EC to obtain g eff A [26][27][28], 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 [29]. In this way, the nuclear shell model and QRPA typically reproduce experimental 2νββ rates, and predict non-measured ones [17,18,[30][31][32].
Reference [39] 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 oddodd 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 [35], and compare the measured T 2ν 1/2 and ξ 2ν 31 values with the predictions from the QRPA and nuclear shell model. Since 0νββ NMEs also show a competition between contributions from lowand 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-µmthick transparent nylon film and is suspended at the center of the KamLAND detector [40,41]. The IB is surrounded by 1 kton of liquid scintillator (LS) which acts as an active shield. The scintillation photons are viewed by 1,879 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/liter of the fluor PPO (2,5-diphenyloxazole), and (2.91 ± 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 ± 0.08)% 136 Xe, (8.96 ± 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. [35] with a total live time of 534.5 days. In order to avoid systematic uncertainties arising from backgrounds, we apply a tightened 2νββ event selection for this work. The fiducial volume for the reconstructed event vertices is defined as a 1-mradius spherical shape at the detector center, which gives a fiducial exposure for this analysis of (126.3 ± 3.9) kgyr in 136 Xe. We perform a likelihood fit to the binned energy spectrum of the selected candidates between 0.8 and 4.8 MeV. The systematic uncertainties on the 2νββ rate are evaluated identically as in Ref. [35] and are summarized in Table I. A detailed energy calibration is essential for the ex-  traction 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 highstatistics 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 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 0νββ and major backgrounds in the Xe-LS, such as 85 Kr, 40 K, 210 Bi, and the 228 Th-208 Pb sub-chain 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 sub-chain 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 [35].
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% 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 non-zero ξ 2ν 31 . The best-fit total 2νββ rate 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. Considering the systematic uncertainties in Ta with energy denominator ∆ = [E j −(E i +E f )/2]/m e . E k is the energy of the nuclear state |J π k 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 [42]. We reproduce M 2ν GT = 0.064 from Ref. [24] with the GCN interaction [19], and also use the alternative MC interaction from Ref. [43], which yields M 2ν GT = 0.024. Both interactions have been used in 0νββ decay studies [11,44]. 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 = q g 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 [24], 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νββ half-life 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 g eff  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 [29]. 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 [26,45].
To illustrate the origin of the differences between the theoretical calculations, Fig. 4 compares the nuclear shell model and QRPA Argonne running sums of M 2ν GT and M 2ν GT −3 [24,39], multiplied by the corresponding (g eff A ) 2 . The sums run over the excitation energy of the spinparity 1 + states in the intermediate nucleus 136 Cs. The theoretical M 2ν GT running sums differ: while the shell model converges at E exc 8 MeV, QRPA terms contribute until E exc 20 MeV. Moreover, at E exc ∼ 10 MeV the accumulated QRPA M 2ν 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 shell-model contributions may be too small due to missing spin-orbit 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,46], 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 high-energy states similar to 2νββ decay is expected [15,47,48]. 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. 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. However, the quenching may not be the same in 2νββ and 0νββ decays, especially in the light of the differences in the two-body [49][50][51] and contact [52] 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 measurement 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. Our analysis reveals that ξ 2ν 31 is sensitive to competing contributions to the NME from low-and highenergy intermediate-states. Since a similar competition is also relevant for 0νββ decay, studies of this observable provide new insights for identifying reliable 0νββ NMEs.
We thank P. Vogel for useful discussions.