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

Abstract

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 ($\hbox {Li}_2^{\;\;100}\hbox {MoO}_4$) 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 (France). From a total exposure of 42.235 kg$\times$day, the half-life of $^{100}$Mo is determined to be $T_{1/2}^{2\nu }=[7.12^{+0.18}_{-0.14}\,\mathrm {(stat.)}\pm 0.10\,\mathrm {(syst.)}]\times 10^{18}$ years. This is the most accurate determination of the $2\nu \beta \beta$ half-life of $^{100}$Mo to date.

1 Introduction

Two-neutrino double-beta ($2\nu \beta \beta$) decay and the related process of the two-neutrino double-electron capture (2$\nu$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\nu \beta \beta$ and to benchmarking the calculations of the beyond the Standard Model process, zero-neutrino double-beta ($0\nu \beta \beta$) decay. The observation of the latter process would demonstrate that the lepton number is not a conserved quantity in nature, and establish the Majorana nature of neutrinos. The search for $0\nu \beta \beta$ is 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_{\beta \beta } = 3034.36(17)$  keV) [4] and the shortest $2\nu \beta \beta$ half-life [3]. Table 1 summarizes the measurements of the $^{100}$Mo  $2\nu \beta \beta$ 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_{1/2}^{2\nu }=[6.81\pm 0.01(\mathrm {stat.})^{+0.38}_{-0.40}(\mathrm {syst.})]$ $\times 10^{18}$ year [5]. The separate foil and detector design limits the scalability of the experiment to large isotope masses; most of the leading $0\nu \beta \beta$ experiments are moving towards the combined detector and isotope design. The most precise previous measurement using the “source = detector” approach is $T_{1/2}^{2\nu }=[7.15\pm 0.37(\mathrm {stat.})\pm 0.66(\mathrm {syst.})]\times 10^{18}$ year [6] using zinc molybdate ($\hbox {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, $\hbox {Li}_2^{\;\;100}\hbox {MoO}_4$ crystals operated as scintillating bolometers and developed as part of the CUPID-Mo program Refs.  [7,8,9,10,11] are used to precisely measure the $^{100}$Mo  $2\nu \beta \beta$ half-life. CUPID-Mo is a demonstrator experiment for CUPID [12, 13], a proposed next-generation bolometric search for $0\nu \beta \beta$ in $^{100}$Mo at the Laboratori Nazionali del Gran Sasso (LNGS, Italy). CUPID will use the infrastructure built for the CUORE $0\nu \beta \beta$ experiment [14], currently in operation at LNGS.

Table 1 Measurements of $2\nu \beta \beta$ decay of $^{100}$Mo to date

2 Experiment

This measurement uses four lithium molybdate crystals enriched in $^{100}$Mo   ($\hbox {Li}_2^{\;\;100}\hbox {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 $\hbox {Li}_2\hbox {MoO}_4$ based scintillating bolometers is described in [7, 8, 25]. Reference [7] contains a conservative estimation of the error in the enrichment (±0.2%) for the two samples produced from the first crystal boule. We have improved this figure for the current study, considering that the $^{100}$Mo enrichment of the molybdenum oxide precursor varied slightly among the different batches of powder used for crystal productions – with typical enrichment uncertainties of the order of ±0.05% – and taking into account 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.

Table 2 $\hbox {Li}_2^{\;\;100}\hbox {MoO}_4$ crystal scintillators used in the experiment

Each $\hbox {Li}_2^{\;\;100}\hbox {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 $\oslash 44\times 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 $\gamma$($\beta$) and $\alpha$ events [28]. This discrimination is key in the analysis that follows for both the estimation and reduction of backgrounds. In fact, light-assisted particle identification allows us to eliminate energy-degraded $\alpha$ particles from surface contamination that leak in the spectral region of interest. It also helps determine the internal contamination of the crystals through $\alpha$ spectroscopy, an important ingredient of a background model. In addition, the detector output is affected by fake instrumental events, probably related to stress releases in materials thermally coupled to the detector (glue or Kapton^® pads) similar to those discussed in Ref. [29]. These events that can be effectively eliminated not only by pulse shape discrimination but also because they are associated with no light signal.

Fig. 1

Left: EDELWEISS-III cryostat modeled in Geant4 as configured for the experiment. Bottom right: zoom on the detector region. Top right: zoom on an individual detector, consisting of a crystal in a copper holder that also supports the light detector, based on a germanium wafer; the details of the copper assembly surrounding the crystals are not shown

The experiment operated in the low-background cryostat of the EDELWEISS-III dark-matter experiment [30], see Figs. 1 and 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 EDELWEISS-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 composition. 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 $\approx$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 $\pm 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 non-thermal variations in the detector baselines. As a result, $\sim$7% and $\sim$12% of physics data were not considered in the present analysis for setups 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\%$.

Fig. 2

EDELWEISS cryogenic facility with partially installed detector modules

The energy scale and energy resolution of the detectors are calibrated using natural radioactivity lines of $^{40}$K and $^{232}$Th, as well as a periodically deployed $^{133}$Ba source [7, 8]. The energy resolution is measured at 356.0 keV ($^{133}$Ba), 1460.8 keV ($^{40}$K) and 2614.5 keV ($^{208}$Tl) resulting in $\sim 3$ keV, $\sim 5$ keV and $\sim 6$ keV full width at half maximum (FWHM), respectively. The energy scale is stable to within $\pm 0.12$% as determined from the variation observed in the periodic $^{133}$Ba calibrations and the physics data ($^{210}$Po $\alpha$ events originating in the crystal bulk). After applying the linear energy calibration based on the position of the $^{208}$Tl 2615 keV line, we observe a modest residual non-linearity in the detector response, manifested as $\pm 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 [31].

The energy spectra of events acquired in setup 1 and setup 2 are shown in Fig. 3. The slight difference between the two spectra, especially at low energies, can be related to some contamination nearby crystal No.1 in setup 1 as this detector was re-assembled using new materials for setup 2. The probability of coincidence between events in the crystals is small due to the detector positions in the setup and a $\sim$2-mm-thick copper shield surrounding each detector. Thus, any coincidences are ignored (i.e. no anti-coincidence cut is applied) in the analysis and in the simulations. A pulse-shape discrimination cut is applied to the signals to find physical events and to reject pileup events; this reduces tails in the energy resolution function. In addition, $\alpha$ decays are eliminated from the spectrum using the light-assisted particle identification cut, which achieves about 9$\sigma$ $\alpha$/$\gamma$ separation [7, 8]. The light-assisted particle identification removes not only fully contained $\alpha$ events from U/Th chains, expected above 4 MeV, but also $\alpha$ decays with degraded energies originating near the crystal surfaces. The rate of energy-degraded $\alpha$ events is estimated to be $0.1-0.2$ counts/year/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\pm 1.2)\%$ and $(96.6\pm 0.7)\%$ for setups 1 and 2, respectively. The exposure-weighted efficiency for the complete data set is $(96.5\pm 0.6)\%$. The selection efficiency estimate was cross-checked using a prominent, but still low intensity, $\gamma$ peak of $^{40}$K resulting in a good agreement at $(94.7\pm 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 $\beta$) to measure the selection efficiency of $(90\pm 10)\%$ at low energies.

Fig. 3

Energy spectra accumulated with the $\hbox {Li}_2^{\;\;100}\hbox {MoO}_4$ scintillating bolometers in the setup 1 (solid blue histogram, exposure 10.308 kg$\times$day) and setup 2 (red dotted histogram, exposure 31.927 kg$\times$day). The $\gamma$-ray energies are listed in keV

3 Results and discussion

3.1 Background model

The most striking feature in the summed energy distribution in Fig. 3 is the continuous spectrum characteristic of $2\nu \beta \beta$ decays, which dominates the data above $\sim 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, which is broader than the lines of similar and higher energies, is consistent with the presence of two sources: an external one, far from the detectors, that contributes with a $\gamma$ ray of 1461 keV only (emitted after the EC decay of $^{40}$K), and an internal one that deposits an additional energy of 3.2 keV corresponding to the K-shell binding energy released after the primary EC decay. The ratios of the other peaks to the continuum indicates that the $\gamma$-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 $\alpha$ region of the energy spectrum.

Based on these observations, we construct a comprehensive background model which includes a combination of “internal” sources (inside the $\hbox {Li}_2^{\;\;100}\hbox {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 $\beta$ events). The backgrounds are simulated using the Geant4 package version 10.p03 (Livermore physics list) [32,33,34] with initial kinematics given by the DECAY0 event generator [35, 36]. The following “external” sources are simulated on the 300 K cryostat vessel indicated in Fig. 1:

• $^{40}$K;

• $^{228}$Ac;

• 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 EDELWEISS setup [30, 37].

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 [30].

The following “internal” sources are simulated inside the crystals:

• $^{40}$K;

• $^{87}$Rb;

• $^{90}$Sr and $^{90}$Y;

• $^{210}$Pb/$^{210}$Bi;

• $2\nu \beta \beta$ decay of $^{100}$Mo to the ground state of $^{100}$Ru;

• $2\nu \beta \beta$ decay of $^{100}$Mo to the first excited state of $^{100}$Ru, $0^+$ at 1130.3 keV. The half-life of this decay is fixed to the value determined by the NEMO-3 Collaboration [38].

The $^{210}$Pb/$^{210}$Bi contribution is determined by the analysis of the $^{210}$Po peaks in the $\alpha$-decay region of the energy spectrum and by taking into account the time elapsed from the growth process for each crystal. 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 $\beta$ 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 [7], and similarity of lithium, potassium and rubidium chemical properties. The presence of $^{90}$Sr-$^{90}$Y in the crystals cannot be excluded.

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 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 $\gamma$-peaks.

It should be stressed that no $\alpha$ decays from the U/Th chain, but few tens–hundreds $\mu$Bq/kg of $^{210}$Po, were observed in the $\hbox {Li}_2^{\;\;100}\hbox {MoO}_4$ crystal scintillators [7, 8], 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.

Table 3 Radioactive contamination of $\hbox {Li}_2^{\;\;100}\hbox {MoO}_4$ crystal scintillators. The limits are quoted at 90% CL

We constrain the background model and the $2\nu \beta \beta$ half-life by performing an extended maximum-likelihood fit [39] to the sum spectrum (the total exposure is 42.235 kg$\times$day, or $3.798(9)\times 10^{23}$ $^{100}$Mo nuclei$\times$ year), binned uniformly with 20 keV bins. We perform a complementary binned least-squares/maximum likelihood fit using PAW/MINUIT software [40, 41]; 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\nu \beta \beta$ 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\nu \beta \beta$ decay, describes the experimental data in the [120–3000] keV range with $\chi ^2=121$ for 126 degrees of freedom.

Fig. 4

Bottom: The energy spectrum accumulated with $\hbox {Li}_2^{\;\;100}\hbox {MoO}_4$ scintillating bolometers (exposure is 42.235 kg$\times$day) and the fit in the energy range [120–3000] keV. The data points represent the data, the solid blue line shows the sum of all components of the fit, the solid red line is the $2\nu \beta \beta$ contributions, and the other components of the fit are described in the legend. Top: fit residuals normalized by the statistical error of the data in each energy bin (“pulls”). The pulls are shown only for the data in the fit range of [120–3000] keV

3.2 Model of the $2\nu \beta \beta$ decay

We simulate $2\nu \beta \beta$ 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\nu \beta \beta$ decay was proposed in [42] for nuclei where the $1^+$ ground state of the intermediate nucleus may dominate the $2\nu \beta \beta$ decay. $^{100}$Mo is one of a few cases where the SSD mechanism is expected to have some merit [43,44,45,46,47]. The data of the NEMO-3 experiment favor the SSD mechanism in $^{100}$Mo [5, 48, 49] and are inconsistent with the HSD hypothesis.

The energy spectra of single electrons and summed two-electron energy spectra for the $^{100}$Mo$\rightarrow ^{100}$Ru $2\nu \beta \beta$ decay using calculations with the SSD and the HSD approximations [47] 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\nu \beta \beta$ decays in the baseline fit to the experimental data, treating the difference between HSD and SSD models as a systematic uncertainty (see Sect. 3.4).

Fig. 5

Single-electron spectra (a) and summed energy spectra of two electrons (b) for the $^{100}$Mo$\rightarrow ^{100}$Ru $2\nu \beta \beta$ decay calculated in the HSD and SSD models. The spectra are normalized to unit area

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\nu \beta \beta$ decays. We perform the fit in the interval [120–3000] keV range using the spectra generated under the SSD and HSD hypotheses. The quality of both fits is acceptable, but the HSD hypothesis returns a larger overall $\chi ^2$ by 12.5 units (the negative log-likelihood is larger for the HSD hypothesis by 8.2 units) .

Since $\sqrt{\Delta \chi ^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 (\mathcal {L_\mathrm {HSD}}/\mathcal {L_\mathrm {SSD}})$ ($\mu =+8.04$), its standard deviation ($\sigma =2.68$), and the probability for the ratio $-\log (\mathcal {L_\mathrm {HSD}}/\mathcal {L_\mathrm {SSD}})$ to fluctuate below zero ($p=0.0014\pm 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\sigma$.

3.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\nu \beta \beta$, the correlations between the background source activities and the $2\nu \beta \beta$ half-life are significant. When fit over the broad energy range, e.g. [120–3000] keV, the best-fit $2\nu \beta \beta$ 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 $\gamma$ 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 $\gamma$ spectra from external sources, and a possible contribution from internal $^{90}$Sr-$^{90}$Y $\beta$-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\nu \beta \beta$ spectrum. 9183 events are found in this range in the 42.235 kg$\times$day of exposure, with 91% attributed to $2\nu \beta \beta$ events. The fit quality is excellent ($\chi ^2=50$ for 61 degrees of freedom) with modest (80%) correlations between the $2\nu \beta \beta$ half-life and the background components. The fit returns $8370^{+162}_{-214}\,(\mathrm {stat.})$ $2\nu \beta \beta$ events in the fit region (extrapolated to the full energy range, the number of $2\nu \beta \beta$ events is $35638^{+693}_{-912}\,(\mathrm {stat.})$). Taking into account the selection efficiency ($0.9646\pm 0.0060$), we find the half-life $T_{1/2}^{2\nu }=[7.12^{+0.18}_{-0.14}\,(\mathrm {stat.})] \times 10^{18}$ year. The statistical 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. The systematic errors are discussed in detail in the following section.

For comparison, the energy interval [120–3000] keV contains 63,717 events; the fit attributes $35405\pm 605$ to $2\nu \beta \beta$ ($99.4\%$ of the $2\nu \beta \beta$ spectrum is contained in the [120–3000] keV interval). We find $T_{1/2}^{2\nu }=[7.13\pm 0.12\,(\mathrm {stat.})\pm 0.20\,(\mathrm {syst.})]\times 10^{18}$ year 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.

3.4 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_{1/2}^{2\nu }$. Uncertainty in the energy scale contributes $0.2\%$. Variation of the bin width (from 10 keV to 30 keV) change $T_{1/2}^{2\nu }$ by up to $0.8\%$. We attribute this variation to the uncertainty in the energy 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 $\hbox {Li}_2^{\;\;100}\hbox {MoO}_4$ crystal scintillators by U/Th is very low. Assuming the activities of the $\beta$ 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 $\le 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 ($\le 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_{1/2}^{2\nu }$ 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 $\gamma$-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 ($\chi ^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_{1/2}^{2\nu }$ by $0.4\%$; we consider this difference to be a conservative upper limit on the systematic error induced by the uncertainty in $2\nu \beta \beta$ spectral shape. The description of the $2\nu \beta \beta$ 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\nu \beta \beta$ half-life, our $2\nu \beta \beta$ decay model does not currently include ${\mathcal {O}}(\alpha )$ and ${\mathcal {O}}(\alpha Z)$ radiative corrections other than the Coulomb (final state) corrections computed in Ref. [47].

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\nu \beta \beta$ decay to the first $0^+$ excited level of $^{100}$Ru. The summary of the systematic uncertainties is given in Table 4.

Table 4 Estimated systematic uncertainties (%)

4 Summary

Adding all systematic contributions in quadrature, the half-life of $^{100}$Mo relative to the $2\nu \beta \beta$ decay to the ground state of $^{100}$Ru is:

$$\begin{aligned} T_{1/2}^{2\nu }=[7.12^{+0.18}_{-0.14}\,(\mathrm {stat.})\pm 0.10\,(\mathrm {syst.})]\times 10^{18}~\mathrm {year}. \end{aligned}$$

This that can be simplified further by summing in quadrature the systematic and statistical errors:

$$\begin{aligned} T_{1/2}^{2\nu }= (7.12^{+0.21}_{-0.17})\times 10^{18}~\mathrm {year}. \end{aligned}$$

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).

Fig. 6

A historical perspective of $T_{1/2}^{2\nu }$ of $^{100}$Mo as a function of the publication date in the experiments: (1) ELEGANT V [15], (2) NEMO-2 [16], (3) segmented Si(Li) detector [18], (4) NEMO-2 reanalyzed [17], (5) Hoover Dam using a time-projection chamber [19], (6) DBA (liquid argon detector) [20], (7) geochemical experiment [21], (8) preliminary result of NEMO-3 [22], (9) low temperature $\hbox {ZnMoO}_4$ bolometers [6], (10) final result of NEMO-3 [5], (11) present study

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\nu \beta \beta$ decay of $^{100}$Mo to the ground state of $^{100}$Ru, assuming the SSD mechanism, can be calculated as $|M^\mathrm {eff}_{2\nu }|=0.184^{+0.002}_{-0.003}$ by using the phase-space factor $4134\times 10^{-21}$ $\hbox {year}^{-1}$ calculated in [47].

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 [38]) is:

$$\begin{aligned} T_{1/2}=(7.05^{+0.21}_{-0.17})\times 10^{18}~\hbox {year}. \end{aligned}$$

In other words, the branching ratios are $99.06(11)\%$ and $0.94(11)\%$ for the $2\nu \beta \beta$ decay of $^{100}$Mo to the ground state and to the first $0^+$ 1130.3 keV excited level of $^{100}$Ru, respectively.

5 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 $\approx 0.2$ kg each operated in the EDELWEISS-III low-background setup at the Modane underground laboratory (France). The $^{100}$Mo half-life value $T_{1/2}^{2\nu }=7.12^{+0.18}_{-0.14}\,(\mathrm {stat.})\pm 0.10\,(\mathrm {syst.})]\times 10^{18}$ year is measured with 42.235 kg$\times$day 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.

Moreover, the half-life value measured with the relative uncertainty of $^{+2.9}_{-2.4}\%$ is among the most precise measurements of any $2\nu \beta \beta$ decay to date. Other leading measurements are of $^{130}$Te by CUORE (2.8% [56]), $^{136}$Xe by EXO-200 (2.8% [57]) and KamLAND-Zen (3.3% [58, 59]), $^{76}$Ge by GERDA (4.9% [60]), $^{116}$Cd by Aurora (5.3% [61]), $^{82}$Se by NEMO-3 (6.4% [62]) and CUPID-0 ($^{+2.2}_{-1.6}\%$ [63]), $^{150}$Nd by NEMO-3 (7.1% [64]), $^{96}$Zr by NEMO-3 (8.9% [65]) and other observations of $2\nu \beta \beta$ decay in $^{48}$Ca, $^{128}$Te, and $^{238}$U ($\approx$10–30%, e.g. see in [3]). The three of four most precise measurements of the $2\nu \beta \beta$ half-life are from the bolometric experiments, demonstrating the power of the technique.

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 $\hbox {Li}_2^{\;\;100}\hbox {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\nu \beta \beta$ measurements.

In agreement with the observation by NEMO-3 [5, 48, 49], we favor the SSD mechanism of the $2\nu \beta \beta$ decay over the HSD mechanism, with the statistical significance of $>3\sigma$. 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 [11] running now in its first phase with 20 enriched $\hbox {Li}_2^{\;\;100}\hbox {MoO}_4$ scintillating bolometers (with mass $\approx 0.2$ kg each).