Precise measurement of $2\nu\beta\beta$ decay of $^{100}$Mo with the CUPID-Mo detection technology

We report the measurement of the two-neutrino double-beta ($2\nu\beta\beta$) decay of $^{100}$Mo to the ground state of $^{100}$Ru using lithium molybdate (\crystal) scintillating bolometers. The detectors were developed for the CUPID-Mo program and operated at the EDELWEISS-III low background facility in the Modane underground laboratory. From a total exposure of $42.235$ kg$\times$d, the half-life of $^{100}$Mo is determined to be $T_{1/2}^{2\nu}=[7.12^{+0.18}_{-0.14}\,\mathrm{(stat.)}\pm0.10\,\mathrm{(syst.)}]\times10^{18}$ years. This is the most accurate determination of the $2\nu\beta\beta$ half-life of $^{100}$Mo to date. We also confirm, with the statistical significance of $>3\sigma$, that the single-state dominance model of the $2\nu\beta\beta$ decay of $^{100}$Mo is favored over the high-state dominance model.


Introduction
Two-neutrino double-beta (2νββ) decay and the related process of the two-neutrino double-electron capture (2νECEC) are allowed second-order processes in the Standard Model theory of electroweak interactions and are the rarest nuclear processes ever observed [1,2,3]. Precise measurements of these processes are critical to understanding the nuclear physics governing 2νββ and to benchmarking the calculations of the beyond the Standard Model process, zero-neutrino double-beta (0νββ) decay. The observation of the latter process would establish the Majorana nature of the neutrino and is therefore the subject of a global experimental effort.
Double-beta decay is observable in nuclei where the single beta decay is forbidden or highly suppressed. Of the candidate isotopes, 100 Mo is characterized by one of the largest decay energies (Q ββ = 3034.36 (17) keV) [4] and the shortest 2νββ half-life [3]. Table 1 summarizes the measurements of the 100 Mo 2νββ half-life to date. Most experiments have used 100 Mo foils coupled with traditional tracking and calorimetry techniques. NEMO-3 presents the most precise measurement to date at T 2ν 1/2 = [6.81 ± 0.01(stat.) +0.38 −0.40 (syst.)] yr [14]. The separate foil and detector design limits the scalability of the experiment to large isotope masses; most of the leading 0νββ experiments are moving towards the combined detector and isotope design. The most precise previous measurement using the "source = detector" approach is T 2ν 1/2 = [7.15 ± 0.37(stat.) ± 0.66(syst.)]×10 18 yr [13] using zinc molybdate (ZnMoO 4 ) crystals operated as scintillating bolometers.
The bolometric technique is now competitive with the foil-based detectors, and offers distinct advantages. In this work, Li 2 100 MoO 4 crystals operated as scintillating bolometers and developed as part of the CUPID-Mo program [15,16,17,18,19] are used to precisely measure the 100 Mo half-life. CUPID-Mo is a demonstrator experiment for CUPID [20,21], a proposed next-generation bolometric search for 0νββ in 100 Mo at the Laboratori Nazionali del Gran Sasso (LNGS). CUPID will use the infrastructure built for the CUORE 0νββ experiment [22], currently in operation at LNGS.

Experiment
This measurement uses four lithium molybdate crystals enriched in 100 Mo (Li 2 100 MoO 4 ) instrumented as low temperature scintillating bolometers. The crystals were produced as part of the LUMINEU project [23]. They were grown using the low-thermal-gradient Czochralski technique starting from the highly purified enriched molybdenum oxide and lithium carbonate [24]. The R&D of large volume Li 2 MoO 4 based scintillating bolometers is described in [16,15,25]. The 100 Mo enrichment of the molybdenum oxide precursor varied slightly between the different crystal productions. The final enrichment fraction in the grown crystals and its uncertainty was estimated taking into account the uncertainties in the original precursor enrichment (±0.05%) and the effect of the crystal growth process. Table 2 summarizes the crystal dimensions, masses, isotopic abundance of 100 Mo in the crystals, and number of 100 Mo nuclei.
Each Li 2 100 MoO 4 crystal is instrumented with a neutron transmutation doped (NTD) germanium temperature sensor [26] and a heavily-doped silicon heater. The latter is used to stabilize the thermal response of the detector [27]. The two devices are glued to the crystal surface and then the crystals are installed in a copper holder and secured by PTFE support clamps. A light detector constructed from a Germanium disc ⊘44 × 0.17 mm instrumented with an NTD sensor is installed above each crystal to detect the scintillation signal from the crystal. The simultaneous detection of heat and light signals provides a powerful discrimination between γ(β) and α events [28]. This discrimination is key in the analysis that follows for both the estimation and reduction of backgrounds.
The experiment operated in the low-background cryostat of the EDELWEISS-III dark-matter experiment [29], see Figs. 1-2. The cryostat is located in the Modane underground laboratory (France) at the depth of 4800 m of water equivalent. The central volume of the EDEL-WEISS-III cryostat is shielded by 20 cm of Pb, the innermost 2 cm is Roman Pb to reduce the 210 Pb background contribution. The experiment was realized in two steps: a single crystal configuration, "setup 1", and a four-crystal configuration, "setup 2". The detector modules and materials in the two setups were slightly different, producing a somewhat different background com-  position. EDELWEISS germanium detectors were run concurrently with this measurement.
In addition to the differences in geometry and materials between the two detector configurations, there was a change in the data acquisition during setup 2. Setup 1 and ≈22% of setup 2 were triggered online, while the remainder of setup 2 was acquired in the streaming data acquisition mode and then triggered offline. The data acquired during instabilities of the cryogenic system were not used for the analysis. If the temperature of the detector holder plate showed variations larger than ±0.1 mK from a chosen value (20.0 mK and 19.2 mK for setup 1, and 17.0 mK for setup 2), the data were discarded. Similarly, we discard periods of large nonthermal variations in the detector baselines. As a result, ∼7% and ∼12% of physics data were not considered in the present analysis for setup 1 and 2, respectively. Table 2 summarizes the live-time for each configuration. The uncertainty in the live-time calculation is estimated from the loss of the periodically injected heater signals. This uncertainty for the online-triggered data is 0.23% and the uncertainty in the stream mode is 0.22%, leading to the exposure-weighted average of 0.22%.
The energy scale and energy resolution of the detectors are calibrated using 40 K, 133 Ba, and 232 Th gamma sources [16,15]. The energy resolution is measured at 356.0 keV ( 133 Ba), 1460.8 keV ( 40 K) and 2614.5 keV ( 208 Tl) resulting in ∼ 3 keV, ∼ 5 keV and ∼ 6 keV full width at half maximum (FWHM), respectively. The energy scale is stable to within ±0.12% as determined from the variation observed in the periodic 133 Ba calibrations and the physics data ( 210 Po α events originating in the crystal bulk). After applying the energy calibration, we observe a modest residual non-linearity in the detector response, manifested as ±5 keV shifts in the position of the known background peaks in the physics data. We correct for these shifts by applying a 2nd-order polynomial correction to the spectra of the reconstructed energies, binned in 1 keV intervals [30].
The energy spectra of events acquired in setup 1 and setup 2 are shown in Fig. 3. Coincidences between the crystals exist but are neglected in the analysis, taking into account a rather small coincidence probability due to the detector positions in the setup and a thick copper shield (minimum 2 mm) surrounding each detector. A pulse-shape discrimination cut is applied to the sig- nals to find physical events and to reject pileup events; this reduces tails in the resolution function. In addition, α decays are eliminated from the spectrum using the light-assisted particle identification cut, which achieves about 9σ α/γ separation [16,15]. The light-assisted particle identification removes not only fully contained α events from U/Th chains, expected above 4 MeV, but also α decays with degraded energies originating near the crystal surfaces. The rate of such events is estimated to be 0.1 − 0.2 counts/yr/kg/keV in the 2.7 − 3.9 MeV region.
The selection efficiency is found to be constant above 500 keV, and is evaluated to be (96.1±1.2)% and (96.6± 0.7)% for setups 1 and 2, respectively. The exposureweighted efficiency for the complete data set is (96.5 ± 0.6)%. The selection efficiency estimate was cross-checked using a prominent, but still low intensity, γ peak of 40 K resulting to a good agreement: (94.7 ± 1.6)%. Below 500 keV, the raw spectrum of triggered events before the light yield and pulse shape selection has a significant contribution from fake instrumental events. We use events identified as 210 Pb decays (a 46.5 keV x-ray and the corresponding β) to measure the selection efficiency of (90 ± 10)% at low energies.

Background model
The most striking feature in the summed energy distribution in Fig. 3 is the continuous spectrum characteristic of 2νββ decays, which dominates the data above ∼ 1 MeV. The most prominent peaks in the spectra can be ascribed to contamination from 40 K and the daughters of the 238 U and 232 Th decay chains. The observed line shape of the 40 K peak is consistent with a subset of the events containing the coincidence of the primary EC decay and subsequent relaxation energy of the atomic shell with the γ-ray event. This indicates two 40 K sources: an external source far from the detectors, and a source internal to the crystals. The ratios of the other peaks to the continuum indicates that the γ-line activity is dominated by the external, far sources that are partially attenuated by the lead and radiation shields of the cryostat. This conclusion is consistent with the limits on the internal crystal contamination in 238 U and 232 Th obtained from the analysis of the α region of the energy spectrum.
Based on these observations, we construct a comprehensive background model which includes a combination of "internal" (inside the Li 2 100 MoO 4 crystals), "external" sources (e.g. detector support structures and the cryogenic vessels), and "nearby" sources (surfaces close to the crystals, where one may expect a contribution from β events). The backgrounds are simulated using the Geant4 package version 10.p03 (Livermore physics list) [31,32,33] with initial kinematics given by the DE-CAY0 event generator [34,35]. The following "external" sources are simulated on the 300 K cryostat vessel indicated in Fig. 1: late 232 Th chain: 212 Pb, 212 Bi and 208 Tl, assumed to be in secular equilibrium; -late 238 U chain: 214 Pb and 214 Bi in secular equilibrium; -137 Cs, which was observed previously in the EDEL-WEISS setup [29,36].
The following "nearby" sources are simulated in the materials near the detectors: -210 Pb/ 210 Bi, assumed to be in secular equilibrium; -208 Tl in the Kapton-based readout connectors, which are known to have measurable levels of contamination [29].
The following "internal" sources are simulated inside the crystals: The 210 Pb/ 210 Bi contribution is determined by the analysis of the 210 Po peaks in the α-decay region of the energy spectrum. The majority of the 210 Pb/ 210 Bi/ 210 Po contamination is attributed to the bulk of the crystals; this is also supported by the shape of the 210 Pb x-ray and β spectra in the vicinity of 46.5 keV. A small contribution from the "nearby" sources (which appears primarily in setup 1) is treated as a systematic uncertainty. "Internal" contamination of 40 K and 87 Rb in the bulk of the crystals are added taking into account the observation of 40 K in some lithium molybdate crystals [15], and similarity of lithium, potassium and rubidium chemical properties. The presence of 90 Sr-90 Y in the crystals cannot be excluded; it is seen with marginal significance in another bolometric experiment [38].
A possibility of the full background reconstruction in a low background experiment is limited by imprecise knowledge of the locations of radioactive contaminations. We build two models with different assumptions about the localization of the background sources. In the default model, we simulate the full geometry of the EDELWEISS cryostat including its payload, and assign 6 all of the "external" contamination to the 300 K vessel. As a systematic check, we also develop a simplified model in which the radioactive backgrounds are placed in copper shields of different thickness around the crystal. This model is tuned to reproduce the energy dependence of the observed intensities of the γ-peaks.
It should be stressed that no α decays from the U/Th chain, but few tens-hundreds µBq/kg of 210 Po, were observed in the Li 2 100 MoO 4 crystal scintillators [16,15], resulting in the very stringent upper limits given in Table 3. Therefore, bulk U/Th radioactivity of the crystals (except for the contribution of 210 Bi) is ignored in the background model, taking into account that the activity of 100 Mo in the crystals is at least three orders of magnitude higher than the possible activity of U/Th daughters. We constrain the background model and the 2νββ half-life by performing an extended maximum-likelihood fit [47] to the sum spectrum (the total exposure is 42.235 kg×d, or 3.798(9) × 10 23 100 Mo nuclei×yr), binned uniformly with 20 keV bins. We perform a complementary binned least-squares/maximum likelihood fit using PAW/MINUIT software [48,49]; the two software packages return consistent results. The background model describes the data very well over a broad energy range [120 − 3000] keV (Fig. 4). In order to assess the sensitivity of the background model and the 2νββ half-life to the underlying assumptions about the background composition, we vary the energy range of the fit in 20 keV steps from 120 to 2000 keV (starting point) to 2300-3000 keV (final point), and find the value of the half-life stable within the expected statistical variations. The model, assuming the single-state dominance mechanism of the 2νββ decay, describes the experimental data in the [120−3000] keV range with χ 2 = 121 for 126 degrees of freedom.

Model of the 2νββ decay
We simulate 2νββ distributions using two assumptions about the decay mechanism: the closure approximation (in other words, high-state dominance, HSD), and the single-state dominance (SSD) hypothesis. The SSD mechanism of 2νββ decay was proposed in [39] for nuclei where the 1 + ground state of intermediate nucleus may dominate the 2νββ decay. 100 Mo is one of a few cases where the SSD mechanism is expected to have some merit [40,41,42,43,44]. The data of the NEMO-3 experiment favor the SSD mechanism in 100 Mo [14,45,46] and are inconsistent with the HSD hypothesis. The energy spectra of single electrons and summed two-electron energy spectra for the 100 Mo→ 100 Ru 2νββ decay using calculations with the SSD and the HSD approximations [44] are shown in Fig. 5. There is a meaningful difference in the single-electron spectra for the HSD and SSD models at low energies, while in the summed energy spectra, measured by bolometric detectors, the difference is substantially smaller. NEMO-3 analysis of the single-electron spectra in 100 Mo rules out the HSD hypothesis with high significance.
We use the SSD model of the 2νββ decays in the baseline fit to the experimental data, treating the difference between HSD and SSD models as a systematic uncertainty (see Section 3.4).
The high statistics of the dataset, excellent resolution, and a high signal-to-background ratio for energies above 1 MeV allow us to test the spectral shape of the 2νββ decays. We perform the fit in the interval [120 − 3000] keV range using the spectra gener- ated under the SSD and HSD hypotheses. The quality of both fits is acceptable, but the HSD hypothesis returns a larger overall χ 2 by 12.5 units (the negative log-likelihood is larger for the HSD hypothesis by 8.2 units) .
Since ∆χ 2 is not, strictly speaking, equal to the significance of discriminating one hypothesis over another [50], we use an ensemble of 10,000 pseudo-experiments to determine the confidence level at which the SSD hypothesis is preferred over the HSD hypothesis. In each pseudo-experiment, we generate the energy distribution of signal and background events from the probability density functions returned by the fit to the [120 − 3000] keV range. The events are generated using the SSD hypothesis, and then two fits using the SSD and HSD hypotheses are performed. From this ensemble, we determine the mean of the expected distribution of the log-likelihood ratio log(L HSD /L SSD ) (µ = +8.04), its standard deviation (σ = 2.68), and the probability for the ratio log(L HSD /L SSD ) to fluctuate above zero (p = 0.0014 ± 0.0004). Similar values are obtained for an ensemble of pseudo-experiments randomly sampled from the energy spectrum observed in the data. We interpret these results as a preference for the SSD hypothesis over HSD with the statistical significance of > 3σ.

Half-life of 100 Mo
The background model described above is sensitive to the exact composition and location of the background sources. Since several possible background sources have broad energy spectra similar to 2νββ, the correlations between the background source activities and the 2νββ half-life are significant. When fit over the broad energy range, e.g. [120 − 3000] keV, the best-fit 2νββ half-life value has a small statistical uncertainty, but a large systematic uncertainty due to the model of background composition and location, as well as the reconstruction efficiency uncertainty at low energies.
For these reasons, we determine the 100 Mo half-life by fit the spectrum in the reduced energy range [1500 − 3000] keV. In this range, only two background contributions are relevant: the late-chain 232 Th decays from external sources, dominated by the 2615 keV 208 Tl γ line and its Compton continuum and the late-chain 238 U decays from external sources, dominated by 214 Bi and its Compton continuum. For completeness, we include a possible contribution from 228 Ac γ spectra from external sources, and a possible contribution from internal 90 Sr-90 Y β-decays. The max-likelihood values of both of those components are consistent with zero. We also split the 208 Tl component into "external" and "nearby" sources. All background components of the fit are restricted to the physical (positive yield) range.
The interval [1500 − 3000] keV contains 23.5% of the 2νββ spectrum. 9183 events are found in this range in the 42.235 kg×d of exposure, with 91% attributed to 2νββ events. The fit quality is excellent (χ 2 = 50 for 61 degrees of freedom) with modest (80%) correlations between the 2νββ half-life and the background components. The fit returns 8370 +162 −214 (stat.) 2νββ events in the fit region (extrapolated to the full energy range, the number of 2νββ events is 35638 +693 −912 (stat.)). Taking into account the selection efficiency (0.9646 ± 0.0060), we find the half-life T 2ν 1/2 = [7.12 +0.18 −0.14 (stat.)] × 10 18 yr. The uncertainties are asymmetric due to the correlations with the background components that are consistent with zero and are restricted to the physical (positive) yield, most notably 90 Y.
For comparison, the energy interval [120−3000] keV contains 63717 events; the fit attributes 35405 ± 605 to 2νββ (99.4% of the 2νββ spectrum is contained in the [120 − 3000] keV interval). We find T 2ν 1/2 = [7.13 ± 0.12 (stat.) ± 0.20 (syst.)] × 10 18 yr for this interval, in excellent agreement to the fit to the more restricted range. The wide energy interval is susceptible to larger systematic uncertainties (discussed below), so we consider this fit as a cross-check. 8

Systematic uncertainties
We vary the underlying assumptions in the default fit over the energy range [1500 − 3000] keV to evaluate the systematic uncertainties. Signal efficiency contributes 0.6% to the systematic error on T 2ν 1/2 . Uncertainty in the energy scale contributes 0.2%. Variation of the bin width (from 10 keV to 30 keV) change T 2ν 1/2 by up to 0.8%. We attribute this variation to the uncertainty in the resolution function applied to the simulated background spectra, and treat the difference as the systematic error.
As it was already mentioned, the internal contamination of the Li 2 100 MoO 4 crystal scintillators by U/Th is very low. Assuming the activities of the β active daughters of 232 Th ( 228 Ac, 212 Pb, 212 Bi, 208 Tl) and 238 U ( 234m Pa, 214 Pb, 214 Bi, 210 Bi) to be equal to the activity limits (see Table 3), the total contribution of the bulk radioactivity is ≤ 0.1% in the region of the fit. The contribution of cosmic muons was estimated on the basis of the measurements with germanium bolometers by the EDELWEISS collaboration [51,52] and the simulations of the muon induced background in germanium detectors taking into account the muon flux as a function of slant depth [53]. A contribution of cosmic-muons background is estimated to be less than 14 counts (≤ 0.15%). We treat these backgrounds as systematic uncertainty (0.2%). In order to further test the sensitivity to the assumptions about the background composition, we repeat the fit after removing the background components consistent with zero activity ( 90 Sr and 228 Ac). As expected, the value of T 2ν 1/2 determined in the [1500 − 3000] keV interval changes very little (0.1%).
We study the effects of the localization of the sources by comparing fits with two independent sets of simulated spectra: one using the complete EDELWEISS geometry and placing all "external" sources on the 300 K vessel, and a simplified detector geometry with location of the sources tuned to reproduce the energy dependence of the observed intensities of the γ-peaks (0.8%). We test the sensitivity to the temporal and spatial variations in the background conditions by splitting the dataset into five independent subsets of similar exposure: setup 1, and 4 separate crystals in setup 2. The five datasets agree within the statistical uncertainties with the half-life determined from the summed spectrum (χ 2 = 2.6 for 4 degrees of freedom). We conclude that there is no evidence for an additional systematic uncertainty arising from this test [54].
The HSD decay model changes T 2ν 1/2 by 0.4%; we consider this difference to be a conservative upper limit on the systematic error induced by the uncertainty in 2νββ spectral shape. The description of the 2νββ energy spectrum can be refined using the improved formalism of the two-neutrino double-beta decay calculations [55]. We should note that like all other measurements of 2νββ half-life, our 2νββ decay model does not currently include O(α) and O(αZ) radiative corrections other than the Coulomb (final state) corrections computed in Ref. [44].
Finally, we account for uncertainties in the number of 100 Mo nuclei, the live-time of the measurement, finite Monte Carlo statistics, and the rate of the 2νββ decay to the first 0 + excited level of 100 Ru. The summary of the systematic uncertainties is given in Table 4.

Summary
Adding all systematic contributions in quadrature, the half-life of 100 Mo relative to the 2νββ decay to the ground state of 100 Ru is: T 2ν 1/2 = [7.12 +0.18 −0.14 (stat.) ± 0.10 (syst.)] × 10 18 yr. This that can be simplified further by summing in quadrature the systematic and statistical errors: 1/2 = (7.12 +0.21 −0.17 ) × 10 18 yr. The half-life value is in an agreement with all the counting experiments after 1995 (a history of 100 Mo half-lives is shown in Fig. 6).
The precision of the present result is higher thanks to the certain advantages of the CUPID-Mo detection technique based on lithium molybdate scintillating bolometers produced from isotopically enriched 100 Mo. The measurement features a high and accurately defined detection efficiency (particularly, because there is no fiducial volume uncertainty), a high energy resolution that allows building an accurate background model, a very low radioactive contamination of the crystal scintillators and of the EDELWEISS-III cryostat. The very low-background conditions, together with utilization of enriched 100 Mo, allowed us to reach a rather high signal to background ratio (approximately 10:1).
An effective nuclear matrix element for 2νββ decay of 100 Mo to the ground state of 100 Ru, assuming the SSD mechanism, can be calculated as |M eff 2ν | = 0.184 +0.002 −0.003 by using the phase-space factor 4134×10 −21 yr −1 calculated in [44]. The effective nuclear matrix element can be written as product M eff 2ν = g 2 A × M 2ν , where g A is axial vector coupling constant, M 2ν is nuclear matrix element. While the value of M 2ν is almost independent on the g A and can be calculated with a reasonable accuracy, the possible range of g A can be quenched from 1.2694 (the free nucleon value) to 0.6-0.8 [56,57,58,59].
Taking into account that 100 Mo nuclei decay by the two modes: to the ground state and to the first 0 + excited level of 100 Ru, the actual half-life of 100 Mo (using the most accurate measurement of the decay of 100 Mo to the first 0 + 1130.3 keV excited level of 100 Ru [37] is: −0.17 ) × 10 18 yr. In other words, the branching ratios are 99.06(11)% and 0.94(11)% for the 2νββ decay of 100 Mo to the ground state and to the first 0 + 1130.3 keV excited level of 100 Ru, respectively.

Conclusions
The two-neutrino double-beta decay of 100 Mo to the ground state of 100 Ru is measured precisely with four 100 Mo-enriched highly radiopure lithium molybdate scintillating bolometers ≈ 0.2 kg each operated in the EDEL-WEISS-III low background setup at the Modane underground laboratory (France). The 100 Mo half-life value T 2ν 1/2 = 7.12 +0.18 −0.14 (stat.) ± 0.10 (syst.)] × 10 18 yr is measured with 42.235 kg×d exposure. The measurement, performed in the energy range [1500 − 3000] keV is statistics-limited, and can be further improved with more data. The result, being in a good agreement with all previous counting experiments after 1995, is the most accurate determination of the 100 Mo half-life.
The high precision of the measurement is achieved thanks to utilization of enriched detectors with an extremely low level of radioactive contamination, operated in the low background environment deep underground. A rather high signal to background ratio in the energy interval of the analysis is reached. The calorimetric approach, together with an excellent energy resolution of the Li 2 100 MoO 4 detectors, ensured a high, clearly defined detection efficiency, and accurate background reconstruction, that are typically the main sources of systematic error in the 2νββ measurements.
In agreement with the observation by NEMO-3 [14,45,46], we favor the SSD mechanism of the 2νββ decay over the HSD mechanism, with the statistical significance of > 3σ. Therefore, we derive the half-life assuming the SSD mechanism of the decay. An effect of the energy spectra shape due to the different mechanisms of the decay is included in the systematic error of the half-life.
The half-life and the spectral shape accuracy are expected to be further improved in the CUPID-Mo experiment [19] running now in its first phase with 20 enriched Li 2 100 MoO 4 scintillating bolometers (with mass ≈ 0.2 kg each).