The background model of the CUPID-Mo $0\nu \beta \beta$ experiment

Abstract

CUPID-Mo, located in the Laboratoire Souterrain de Modane (France), was a demonstrator for the next generation $0\nu \beta \beta$ decay experiment, CUPID. It consisted of an array of 20 enriched Li$_{2}$ $^{100}$MoO$_4$ bolometers and 20 Ge light detectors and has demonstrated that the technology of scintillating bolometers with particle identification capabilities is mature. Furthermore, CUPID-Mo can inform and validate the background prediction for CUPID. In this paper, we present a detailed model of the CUPID-Mo backgrounds. This model is able to describe well the features of the experimental data and enables studies of the $2\nu \beta \beta$ decay and other processes with high precision. We also measure the radio-purity of the Li$_{2}$ $^{100}$MoO$_4$ crystals which are found to be sufficient for the CUPID goals. Finally, we also obtain a background index in the region of interest of 3.7 $^{+0.9}_{-0.8}$ (stat)$^{+1.5}_{-0.7}$ (syst) $\times ~10 ^{-3}$ counts/$\Delta E_{\text {FWHM}}/\text {mol}_{\text {iso}}/\text {year},$ the lowest in a bolometric $0\nu \beta \beta$ decay experiment.

1 Introduction

Neutrinoless double beta decay $(0\nu \beta \beta )$ is a hypothetical nuclear transition that would occur if the neutrino is its own antiparticle, or a Majorana particle. It consists in the transformation of an even-even nucleus into a lighter isobar containing two more protons accompanied by the emission of two electrons and no other particles, with a change of the total lepton number by two units. Thus, the $0\nu \beta \beta$ signal is a peak in the summed electron energy spectrum positioned at the Q$_{\beta \beta }$ (the energy difference between parent and daughter nuclei) of the transition. The detection of this “matter-creating” process would represent the observation of a new phenomenon beyond the Standard Model [1]. Current best limits for $0\nu \beta \beta$ half-life are of the order of 10$^{24}$–10$^{26}$ year [2,3,4,5,6,7,8].

The Standard Model process, two-neutrino double beta decay, $2\nu \beta \beta$, includes also the emission of two $\bar{\nu }_e$ and conserves lepton number. Unlike $0\nu \beta \beta$ decay, $2\nu \beta \beta$ has a continuous energy spectrum and has been observed in more than ten nuclei with half-lives in the range of 10$^{18}$–10$^{24}$ year [9].

One of the largest challenges in $0\nu \beta \beta$ decay experiments is the control of the radioactive background, that may produce events in the signal energy region. These could mimic the very rare $0\nu \beta \beta$ signal reducing the experimental sensitivity.

During the last 10 years the scintillating bolometer technology has proved that bolometers based on lithium molybdate (Li$_2$MoO$_4$), are very promising detectors for next generation $0\nu \beta \beta$ searches [10, 11]. Scintillating bolometers were developed to reduce the background observed in the current leading $0\nu \beta \beta$ bolometric experiment, CUORE [12]. In CUORE, the background in the region of interest is dominated by surface $\alpha$’s emitted from the copper structure holding the detectors [13]. The array of 988 TeO$_2$ bolometers, installed at the Laboratori Nazionali del Gran Sasso, Italy, has observed a background in the $^{130}$Te region of interest $({\mathcal {Q}}_{\beta \beta }$ = 2527 keV) of $(1.49\pm 0.04 )\times 10^{-2}~\text {counts}/\text {keV}/\text {kg}/\text {year}$ [4, 14, 15]. The next generation experiment CUPID (Cuore Upgrade with Particle IDentification) will drastically reduce the background thanks to the simultaneous readout of heat and light signals. The capability to discriminate $\beta /\gamma$ from $\alpha$ particles with scintillating bolometers relies on the fact that the light emitted in the Li$_{2}$ $^{100}$MoO$_4$ by $\alpha$ particles is about a factor 5 smaller compared to the light emitted by $\beta /\gamma$’s of the same energy [10, 16].

In addition to the particle discrimination, the CUPID strategy to reduce the background relies on the radiopurity of the scintillating crystals and the minimisation of the passive materials [17]. Another key point is that the ${\mathcal {Q}}_{\beta \beta }$ energy value of $^{100}$Mo (3034 keV) is higher than 2615 keV, implying a signal located above the majority of $\gamma$ lines from natural radioactivity.

The CUPID-Mo experiment [11], located in the Laboratoire Souterrain de Modane (LSM) in France, under an overburden of 4800 m water equivalent, was built as a demonstrator experiment for CUPID. It consisted of 20 Li$_{2}$ $^{100}$MoO$_4$ (LMO) scintillating bolometers and 20 Ge light detectors (LDs) for a simultaneous read-out of heat and light. One of the aims of CUPID-Mo was to validate the background predictions for CUPID, in particular the LMO crystal radiopurity and residual $\alpha$ background. LMO radiopurity for the U and Th chains of less than 10 and 3 $\upmu$Bq/kg respectively are needed to meet the CUPID goals [17], this can be validated on a mid-scale with CUPID-Mo.

In this paper we present the background model which describes the background sources in the CUPID-Mo experiment. This model is based on fitting the CUPID-Mo data to detailed Monte Carlo simulations. We show that the residual $\alpha$ background contribution and the radiopurity of the LMO crystals are sufficient to meet the CUPID background goal.

CUPID-Mo was also an important experiment in its own right. In particular, it has set the world-leading limits on the half-life of $0\nu \beta \beta$ decay of $^{100}$Mo to both ground and excited states [6, 11, 18]. The detailed study of the experimental backgrounds in the CUPID-Mo experiment enables a high precision measurement of the $2\nu \beta \beta$ decay rate and allows to disentangle between the Single State Dominance (SSD) or High State Dominance (HSD) mechanisms [19, 20] of the $2\nu \beta \beta$ decay process in $^{100}$Mo. It also provides the basis to study new physics processes outside the Standard Model, which could distort the spectral shape of the $2\nu \beta \beta$ spectrum, such as $0\nu \beta \beta$ decay with emission of Majoron(s), $2\nu \beta \beta$ decay with emission of Bosonic neutrinos, Lorentz invariance violation or sterile neutrinos [21,22,23,24,25,26,27,28,29,30].

2 The CUPID-Mo experiment

CUPID-Mo was installed in the EDELWEISS cryogenic set-up [31] at LSM. The experiment was in operation between March 2019 and June 2020.

Fig. 1

Left: An individual CUPID-Mo bolometer showing the transparent Li$_{2}$ $^{100}$MoO$_4$ crystal, the copper holder, the NTD-Ge thermometer and the Teflon clamps. Right: View of the opposite side of the detector, showing the black light detector, fabricated from Ge wafers [6]

Fig. 2

Left: The CUPID-Mo experiment in the EDELWEISS cryostat. The five towers on the front contain the CUPID-Mo detectors and the EDELWEISS detectors can be seen behind. Right: GEANT4 rendering of the detector chamber in the Monte Carlo simulation geometry. Inside the five towers are placed the LMO crystals, the light detectors, the clamps and the reflective foils (not seen). The readout cables and the structure supporting the towers are indicated

CUPID-Mo used $^{100}$Mo-enriched LMO crystals, where the $^{100}$Mo, the double beta isotope, has been enriched at $\sim$ 97%. The basic detection modules are the crystals coupled to thermal sensors, consisting of a Neutron Transmutation Doped Ge thermistor, NTD. The top and bottom of the crystals are facing light detectors fabricated from Ge wafers, also instrumented with NTDs to readout the scintillation light signal. The crystals are housed in cylindrical copper holders and supported by PTFE pieces, as shown in Fig. 1. A reflective foil (3 M Vikuiti) is installed around the crystals, inside the copper holder, to increase the light collection. The average weight of the CUPID-Mo crystals is 210 g and the total mass is 4.158 kg corresponding to 2.264 kg of $^{100}$Mo.

The array of 20 bolometers is arranged in five towers with four modules each, as shown in Fig. 2. Each tower is suspended by stainless steel springs to mitigate the vibrational noise of the set-up. The signal from the CUPID-Mo detectors are readout with NOMEX cables, copper and constantan wires on Kapton pads. Situated in the same cryostat, the EDELWEISS detectors, visible behind the five CUPID-Mo towers in Fig. 2, are equipped with Kapton pads and MillMax connectors. The detector chamber consists of four copper plates made of NOSV grade copper¹ to support the bolometers, and is able to accommodate 12 detector towers.

The cryostat involves five thermal copper screens, typically referred to as the 10 mK, 1 K, 50 K, 100 K and 300 K stages respectively. The cryostat screens are made of NOSV and CUC2² grade copper. An internal polyethylene (PE) shield, to shield against neutrons produced in the set-up components by $(\alpha$, n) reactions or induced by muons [32], is mounted between the detectors and the internal lead shield, and has a temperature of $\sim 1$ K. An internal lead shield of 14 cm Roman lead [33] is installed inside the cryostat at 1 K, between the detector chamber and the dilution unit (see Fig. 3). Its main purpose is to shield the detectors from radioactive background of the warm electronics, the cold electronics and the connectors and cables at the 1 K stage.

The external shielding closest to the cryostat consists of 20 cm thick lead, with the innermost 2 cm made of Roman lead. The empty space between the lead shield and the outermost thermal screen of the cryostat is flushed with radon depleted air from a radon trapping facility. The average radon level in the air supplied by the facility is 20 mBq/m$^3$ [34]. Following the external lead shield, a 50 cm thick polyethylene shield is used to moderate the radiogenic neutron flux. A plastic scintillator based active muon veto system surrounds the whole experiment for muon tagging [35] (see Fig. 4).

2.1 Performances

CUPID-Mo has shown excellent detector performances in terms of energy resolution $(7.4 \pm 0.4)$ keV FWHM at 3034 keV [6] and $\alpha$ particle rejection > 99.9% [36], demonstrating that the CUPID requirements are within reach. Further details on the CUPID-Mo set-up and performances are given in [36].

Fig. 3

GEANT4 rendering of the CUPID-Mo Monte Carlo simulation geometry, showing the cryogenic set-up

Fig. 4

Visualisation of the EDELWEISS cryostat and shielding as implemented in our MC simulations, we show the cryostat surrounded by the lead shield, the external polyethylene shielding and the muon veto panels. The muon panels are free to move to give a full geometric coverage

3 Experimental data

The aim of our data processing is to convert the raw data stream into three calibrated energy spectra: $\beta /\gamma$ like events with energy deposits in a single crystal $({\mathcal {M}}_{1,\beta /\gamma })$, events with energy deposits in two crystals $({\mathcal {M}}_2)$ and of $\alpha$-like events $({\mathcal {M}}_{1,\alpha })$. These spectra will then be used in a simultaneous fit to extract radioactive contamination values and describe the observed spectra. The algorithms used for the data processing are described in detail in [6] but we will give a summary of the most important steps in the following. We also estimate the detector response parameters (energy resolution, energy bias, efficiencies, light yield) which are needed for post-processing the Monte Carlo spectral shapes.

3.1 Data taking

In this paper, we use the same dataset as in [6] with an exposure of 2.71 kg $\times$ year of LMO corresponding to 1.47 kg $\times$ year of $^{100}$Mo. Our data is acquired as a continuous time-stream and digitized at 500 Hz by the EDELWEISS DAQ [36] and stored at both CC-IN2P3 (France) and NERSC (USA) for offline analysis. We acquire runs, periods of around 10–100 h of stable data taking, of both physics and calibration data, where a calibration source was placed in the vicinity of the experiment. We use regular calibrations with a $^{232}$Th/$^{238}$U source to calibrate the LMO detectors and a high activity $^{60}$Co source, which generates $^{100}$Mo X-rays in the detectors, to calibrate the LDs. We divide the data into twelve periods of $\sim 1$ month of stable data taking, called datasets. We discard three short periods of data $(\approx$ 1 week each) due to the low statistics causing an inability to accurately calibrate this data.

3.2 Data processing

We process our data using the C++ softwares Apollo and Diana [37, 38], first developed for the CUORE experiment and also used by CUORICINO, CUPID-0 and CUPID [39].

A complete description of the data processing can be found in [6]. We identify physics events using an optimal trigger, also used for previous CUPID-Mo analysis. This triggering is used for both the LDs and the LMO bolometers. We then store a 3 s waveform for both LDs and the LMO channels for each triggered events. For each LMO we associate (up to) two light detectors, called side-channels, which correspond to the LDs facing this LMO detector. These are numbered S1/S2 where S1 is the LD with the better detector performance (lower noise and higher detector light yield).

We estimate the amplitude of peaks using an optimal filter, which maximises the signal to noise ratio based on inputs of the known signal shape and spectral noise power density. This is done for all LMO events and also the corresponding LD events on the side-channels.

Next we correct for thermal gain changes and calibrate our data using dedicated $^{232}$Th/$^{238}$U calibration measurements. This calibration is accurate to around $<1$ keV [6] which is sufficient for the binned background model fits. The LDs are calibrated using the dedicated $^{60}$Co calibrations which produces $\sim 17$ keV Mo X-rays.

3.3 Multiplicity

We define coincidences between physics events, where multiple detectors are triggered simultaneously. This provides useful information since events of $0\nu \beta \beta$ decay or $2\nu \beta \beta$ decay to the ground state are very likely $(>75\%$ probability [11]) to deposit energy in only one crystal. However, background events in particular from $\gamma$’s are likely to deposit energy in multiple crystals simultaneously, for example due to Compton scattering in one crystal, or multiple $\gamma$’s from the same decay. We estimate the multiplicity of an event as the number of pulses in different LMO detectors above our analysis energy threshold (set at 40 keV) within a $\pm 10$ ms time window.

3.4 Data selection

Several cuts are used to remove non-physical events (for example noise spikes and cross-talk) or coincidences of two or more pulses generated by events very close in time within the same crystal, called pile-up events. We require that there is only one trigger in the 3 s LMO waveform. We then define a pulse shape discrimination (PSD) cut, described in detail in [6], using a principal component analysis method (PCA). We also define a cut on the pulse rise time and optimal filter based PSD variables³ which help to cut pile-up like events. Details of the choice of the selection cuts is given in [6].

3.5 Particle identification

Since CUPID-Mo is a dual readout experiment we can discriminate $\alpha$ from $\beta /\gamma$ particles. The use of light detectors also allows us to remove background events in which a particle deposits energy on our LDs. We select $\beta /\gamma$ candidate events using the LD signal as following. We normalise the measured LD signals by defining the variable n, as the difference between the measured LD energy, $E_{LD}$, and the mean expected $\beta /\gamma$ LD energy L, normalized by the light band width $(\sigma )$. We compute n for ${\mathcal {M}}_1$ events. As each LD has different characteristics, the calculation is done for each channel (crystal) c, each dataset d and for both side channels s, i.e.:

$$\begin{aligned} n_{c,s,d}= \frac{E_{LD}-L_{c,d,s}(E)}{\sigma _{c,d,s}(E)}, \end{aligned}$$

(1)

with E the measured LMO energy. The parameter n$_{c,s,d}$ has a distribution expected to be centered at zero for $\beta /\gamma$’s and at a value different from zero for $\alpha$ particles. For details on the determination of the mean expected LD energies and its uncertainty, see [6].

For events with two LDs we expect the 2D distribution of $n_{c,1}$ against $n_{c,2}$ to be a bivariate Gaussian. As we observe no clear correlation we place a radial cut on the variable:

$$\begin{aligned} D = \sqrt{n_{c,1}^2+n_{c,2}^2}. \end{aligned}$$

(2)

If only one LD is available the cut is instead placed just on this n$_{c,s,d}$. We chose a cut of $D<4$ to select $\beta /\gamma$ events and call this data spectrum ${\mathcal {M}}_{1,\beta /\gamma }$. We also construct a spectrum of ${\mathcal {M}}_{1,\alpha }$ events comprised of high energy ${\mathcal {M}}_1$ events, $E> 3$ MeV, with no light selection applied. This data comprises almost entirely $\alpha$ particles. The same events are obtained with a selection cut $D>4$, thus, for simplification we have chosen only the energy cut to select $\alpha$ events.

Unlike most other analysis of scintillating bolometers we also develop a light selection cut for ${\mathcal {M}}$ $_2$ events as described in detail in [18]. For a ${\mathcal {M}}$ $_2$ event the scintillation light recorded can be the sum of that from the crystal above and below a given LD. We use the modeling described in [6] to compute the expected light detector energy for each physics event accounting for multiple contributions to the light yield. From this we can define the normalised LD energy for each pulse in a ${\mathcal {M}}>1$ event. We require that each normalised LD energy (for each channel and side channel) is between $-10$ and 10 $\sigma$, for all but one LD. In this particular LD we observe an accidental contamination of $^{60}$Co. Therefore we generally observe $\gamma$ events in the LMO and $\beta$ events in the LD with very large energy compared to scintillation light. For the two LMOs adjacent to this LD we make a cut of $-10$ to 3 $\sigma$, to take into account the energy directly deposited in the light detector. For more details see [18]. In addition, to further suppress events from this localized $^{60}$Co source we make a global LD anti-coincidence cut to remove the $\gamma$ background originating from this LD. We remove any events (on non-adjacent LMOs) with a trigger on this LD with energy $>2$ keV within a 5 ms window.

3.6 Muon veto anti-coincidence

Despite the large rock overburden at LSM, which suppresses most muon events, they still form a possible background source. The EDELWEISS cryostat has a muon-veto system to remove these events, as shown in Fig. 4. We remove events, in each of the ${\mathcal {M}}$ $_{1,\beta /\gamma }$, ${\mathcal {M}}$ $_2$ and ${\mathcal {M}}$ $_{1,\alpha }$ spectra, with a trigger in the veto system within a 5 ms window. With 98% geometric coverage and the operation voltage adjusted for the aging of the scintillator we expect an ${\mathcal {O}}$(90%) tagging efficiency of muons with a minimal impact on the $\beta /\gamma$ acceptance [35]. Since this background was already subdominant and is strongly suppressed by the veto cut we do not include muons in our background model.

3.7 Delayed coincidences

Radioactivity from the $^{232}$Th and $^{238}$U decay chains in the LMO crystals could be a significant background in our data. Similar to other analyses of scintillating bolometers [5, 40], we can exploit the time correlation of these decay chain events to reduce our experimental backgrounds. In particular, we veto events from the lower part of both chains where there are backgrounds from $^{214}$Bi $(^{238}$U chain) and $^{208}$Tl $(^{232}$Th chain). For $^{208}$Tl we veto events in $10\times T_{1/2}$ (1830 s) following a suspected $^{212}$Bi $\alpha$-decay. This time window contains $>99.9\%$ of the $^{208}$Tl decays.

The very low CUPID-Mo radioactivity also enables a novel delayed coincidence cut removing $^{214}$Bi candidate events. The $^{222}$Rn decay chain proceeds as follows:

(3)

We can therefore tag the event based on the $^{222}$Rn or $^{218}$Po $\alpha$ events and a fairly long dead time. We use energy cuts of 5985–6145 keV for $^{218}$Po and 5460–5620 keV for $^{222}$Rn to tag $\alpha$ candidates. For either type of $\alpha$ candidate events we then veto events within the same crystal within a time window containing 99% of events which is evaluated with MC sampling as 13860 s for $^{222}$Rn and 13620 s for $^{218}$Po. The two possible cuts on $^{222}$Rn or $^{218}$Po improve the rejection power for surface backgrounds. This cut has a small inefficiency (see Sect. 3.9), despite the long veto time.

3.8 Data spectra

Based on these cuts we construct the three data spectra used in our analysis:

• ${\mathcal {M}}$ $_{1,\beta /\gamma }$: Events in one detector identified as $\beta /\gamma$,

• ${\mathcal {M}}_2$: Events in coincidence between 2 crystals, the two energies deposited in each crystal are summed,

• ${\mathcal {M}}$ $_{1,\alpha }$: Events in one detector with alpha energy scale $(>3$ MeV).

Because of the relatively fast half-life of $2\nu \beta \beta$ in $^{100}$Mo $(\sim 7 \times 10^{18}$ year) and extremely low levels of contamination, relatively few peaks are observed in the ${\mathcal {M}}_{1,\beta /\gamma }$ spectrum, where the spectrum of $2\nu \beta \beta$ decays of $^{100}$Mo is the dominant feature. The secondary datasets, ${\mathcal {M}}$ $_2$ and ${\mathcal {M}}_{1,\alpha }$ however contain a lower fraction of $2\nu \beta \beta$ events and therefore provide useful information to determine the location of radioactive contaminations. The experimental spectra after all cuts are shown in Fig. 5.

Fig. 5

CUPID-Mo experimental data. Left: ${\mathcal {M}}$ $_{1,\beta /\gamma }$: events in one detector identified as $\gamma /\beta$. ${\mathcal {M}}$ $_2$: events in coincidence between 2 crystals, the two energies deposited in the crystals are summed. Right: ${\mathcal {M}}$ $_{1,\alpha }$: events in one detector with $\alpha$ energy scale

3.9 Data selection efficiencies

We evaluate the efficiency of our cuts and correct the MC simulations by these values. In particular we use events in $\gamma$ peaks from ${\mathcal {M}}$ $_2$ and ${\mathcal {M}}$ $_{1,\beta /\gamma }$ spectra to evaluate the efficiency of the PSD (Sect. 3.4), light yield and rise time cuts (Sect. 3.5). We do not observe that the cuts have any energy dependence in the range of the utilised $\gamma$ peaks (236–2615 keV). For cuts where the inefficiency can be considered as a dead time, the multiplicity, muon veto, delayed coincidence and LD anti-coincidence, we evaluate the efficiency using the $^{210}$Po peak. We evaluate the pile-up efficiency, the probability a pulse will be superimposed with another in a 3 s window, using random noise triggers. More details on each of these calculations can be found in [6], and the results are summarised in Table 1.

Table 1 Efficiencies for the cuts used on CUPID-Mo data. The PSD and Light Distance cut efficiencies are evaluated using $\gamma$ peaks, and show no energy dependence in the range of the fit

3.10 Energy scale and resolution

We use the observed $\gamma$ peaks in both background and calibration data to predict the energy linearity and resolution. Each LMO detector in each dataset has a distinct energy resolution. As in [6, 11] we perform a fit of the 2615 keV peak in calibration data to extract the resolution of each detector-dataset pair. We use these resolutions to build a function including a common scale factor R(E) which will be determined for the peaks in background data. For our Monte Carlo simulations for each event we sample from a Gaussian with mean E and width $R\times \sigma _{c,d}$, where $\sigma _{c,d}$ is the energy resolution in channel c and dataset d. This energy calibration is discussed in detail in [6].

3.11 Features of data spectra

We observe in Fig. 5 that the spectrum of $2\nu \beta \beta$ decays of $^{100}$Mo dominates the ${\mathcal {M}}$ $_{1,\beta /\gamma }$ data, whereas the ${\mathcal {M}}$ $_2$ spectrum has significant contributions from natural radioactivity, shown by prominent $\gamma$ peaks. These consist predominantly of decays from the $^{238}$U and $^{232}$Th decay chains, however we also observe contributions from $^{40}$K and cosmogenic activation products $^{60}$Co and $^{57}$Co. We also observe a short lived peak of $^{99}$Mo, present for $\sim 1$ dataset, from neutron activation after a calibration with an AmBe neutron source.

The spectrum ${\mathcal {M}}$ $_{1,\alpha }$ is dominated by $\alpha$ decays from components very close to the detectors. As shown in Fig. 5, in our data we observe a large contribution of $^{210}$Po, $E_{\alpha }=5303$ keV, with both a large Q-value and $\alpha$-energy peak and peaks from several other nuclides in the U/Th chains.

During $\alpha$ decay the energy released is shared between the $\alpha$-particle and recoiling nucleus (NR), with energy ${\mathcal {O}}$(100 keV). In LMO crystals the range of $\alpha$ particles is about 10 $\upmu$m and a few nm for nuclear recoils. Therefore we expect to observe a peak at the Q-value of the decay for a LMO bulk event. For surface activity the energy spectrum depends on the implantation depth. For shallow contribution ${\mathcal {O}}$(nm) in the crystal the $\alpha$ or recoil could escape, or both could be contained in the crystal. We therefore expect peaks at the NR energy, at the $\alpha$ energy, and at the Q-value, with a relatively low flat continuum from partial contained $\alpha$’s or NR. The ratio of these peaks depends on the depth of radioactive contamination. For a deeper contribution ${\mathcal {O}}(\upmu$m) the NR is almost always contained but the $\alpha$’s can still escape after depositing some of its energy, giving rise to a continuum extending from low energies up to the Q-value.

Similarly, for materials facing the crystals we expect a dependence on the implantation depth: at shallow depths the spectrum will be characterised by peaks at the $\alpha$-energy and NR energy, for a deep contribution this will become a flat spectrum from low energy up to the alpha energy. We note from Fig. 5 that we generally do not observe clear $\alpha$-energy peaks in our data. However due to the limited statistics the data is still compatible with a full surface contamination. The lack of clear $\alpha$-energy peaks creates a challenge for assessing the surface contamination.

4 Background sources

The background in our experiment is expected mainly from the natural radioactivity in the whole experimental setup, including the detectors. Other contributions from muons, neutrons and environmental gammas are expected to be subdominant, as explained in Sect. 4.1. To minimize the background, all the materials used to build the experiment have been carefully selected in terms of radiopurity. To this end, the daughters of $^{238}$U and $^{232}$Th decay chains, $^{40}$K, and cosmogenic radionuclides have been measured by High Purity Ge $\gamma$-ray spectroscopy and ICPMS (Inductively Coupled Plasma Mass Spectrometry). The CUPID-Mo materials were chosen to minimize the $^{226}$Ra and $^{228}$Th contaminations, as the most critical radioactive backgrounds in the 3 MeV region relevant to $0\nu \beta \beta$ decay searches arise from $^{214}$Bi and $^{208}$Tl decays. Table 2 reports the radioactivity in the CUPID-Mo detector components resulting from CUPID-Mo and CUORE measurement campaigns [15]. The materials which are directly facing the crystals (all but the springs from Table 2) are referred to as close components in the following. The material choice in the EDELWEISS cryostat was done to minimize the contaminants at lower energies, ${\mathcal {O}}(100$ keV), which is the region of interest in dark matter searches. Table 3 shows the radioactivity in the EDELWEISS cryostat materials.⁴ The NOSV copper is used for the CUPID-Mo detector holders, all the copper parts in the detector chamber and the cryostat screens (with the exception of the 1 K screen).

We identify the most significant contributions to our experimental background using the screening measurements and the analysis of experimental data from Sect. 3.11.

Table 2 Radioactivity values of the components of the CUPID-Mo detectors. All measurements have been made by ICPMS, with the exception of the Springs, measured by $\gamma$-ray spectroscopy and the surface contamination of the Vikuiti^™ reflective foil, measured with the BiPo-3 detector [41]

Table 3 Radioactivity of the components in the EDELWEISS setup. All measurements have been made by HPGe $\gamma$-spectroscopy. The MillMax connectors have also been measured by ICPMS

We can broadly categorise our background sources into four groups:

• Close source: Radioactivity in the LMO crystal, reflective foils, LDs, PTFE clamps and NTDs, directly facing the crystals;

• 10 mK source: Sources of activity in the 10 mK stage of the cryostat but not directly facing the LMO crystals (springs, cables, connectors, copper plates for bolometer support), as shown in Fig. 2;

• Infrastructure source: The copper cryostat screens and the internal shieldings, see Fig. 3;

• External: Activity originating from outside the 300 K Cu shield.

4.1 Other contributions – muons, neutrons and environmental gammas

The muon flux at the LSM is 5 muons/m$^2$/day [35]. Muons would generally deposit energy in multiple detectors and be strongly suppressed by anti-coincidence with the muon veto detector (see Sect. 3.6), therefore we do not include them in the background model.

Neutrons may induce background in the $0\nu \beta \beta$ region of interest (ROI) if they are captured in the materials of the setup, producing high energy gammas. The thermal neutron flux in LSM has been measured as $(3.6\pm 0.05~(\text {stat.})\pm 0.27~(\text {syst.}))\times 10^{-6}$ neutrons/s/cm$^2$ [43] and the ambient neutron flux (fast plus thermal) has been estimated $\sim 10^{-5}$ neutrons/s/cm$^2$ in [43, 44]. Previous work [45] showed that 48 cm of polyethylene reduces the neutron flux by a factor $2 \times 10^6$. Taking into account the surface of the CUPID-Mo detectors, we get that the neutron flux expected is less than 1 neutron/year. Thus, ambient neutrons are not taken into account among our background sources.

The gamma flux at LSM has been measured with a portable Ge detector at several locations in the laboratory. At the place where the EDELWEISS set-up is installed, the flux of 2.6 MeV photons was measured as $5.1 \pm 0.2 \ (\text {stat.})\times 10^{-2}$ $\gamma$/s/cm$^2$ [46]. Considering that 20 cm of lead reduce the flux by about a factor $10^4$, then the contribution of environmental gammas may not be negligible. We expect about 6 photons of 2.6 MeV on the detectors surface during the course of all data taking. We take them into account by generating decays at the level of the outermost cryogenic thermal shield, as the spectral shapes measured in the detector from a source generated outside the external lead and outside the outermost cryogenic thermal differ slightly only below 500 keV.

5 Monte Carlo simulations

The Monte Carlo simulation is developed in GEANT4 and implemented with version 10.04 [47]. The MC simulation program developed by the EDELWEISS collaboration [32], has been adapted to include the CUPID-Mo detectors and to include the features described below.

We generated $2\nu \beta \beta$ decay events, with energies sampled from the theoretical two-dimensional single electron energy spectrum from [48, 49]. We consider separately both the HSD and SSD mechanisms.

The radioactive decays in the components of the experiment are generated using both Decay0 [50] and GEANT4. For decay chains in close sources we use the GEANT4 class G4RadioactiveDecay. This allows to generate sub-chains, for example $^{226}$Ra to $^{210}$Pb. We store the final position of the nuclear recoil, and use it as the initial condition for the next decay, along with the time difference. This allows for the accurate consideration of pile-up from events out of the same decay chain and of the delayed coincidence cuts (see Sect. 3.7). From the simulations we then store the energy deposited in both the LMO and LDs. For the $^{232}$Th decay chain, we generate separately $^{232}$Th, $^{228}$Ra to $^{228}$Th and $^{228}$Th to $^{208}$Pb, since $^{228}$Ra and $^{228}$Th have long half-life and so secular equilibrium cannot be assumed. Similarly for the $^{238}$U chain we generate separately $^{238}$U to $^{234}$U, $^{234}$U, $^{230}$Th, $^{226}$Ra to $^{210}$Pb and $^{210}$Pb to $^{206}$Pb.

We use Decay0 for most external sources, not directly facing the crystals, where pile-up events in the same crystal from subsequent decays in a chain are unlikely, and delayed coincidence cuts through the tagging of $\alpha$ events is impossible. For the $^{238}$U decay chain we generate the $\beta$ emitters $^{214}$Pb and $^{214}$Bi. Since they are in secular equilibrium, we combine their spectra to reduce the number of components in the background model fit. We also generate in some components $^{210}$Bi which is not assumed to be in equilibrium. For the $^{232}$Th decay chain we generate $^{212}$Pb, $^{212}$Bi and $^{208}$Tl out of the $^{228}$Th sub-chain and combine them into one spectrum. We also generate $^{228}$Ac which is not assumed to be in equilibrium.

In addition to the $^{238}$U and $^{232}$Th chains, we simulate $^{40}$K contamination in the springs and the outermost cryogenic thermal shield. We have also considered $^{60}$Co from cosmogenic activation in copper as well as $^{87}$Rb and $^{90}$Sr+$^{90}$Y in the crystals. A $^{99}$Mo contribution due to neutron activation in the first days of data taking is also simulated in the crystals.

The decays are generated in the bulk of the components and also the surface for close sources, where surface contaminants can produce a distinct energy spectrum compared to bulk contamination. Surface contaminations are modelled with an exponential density profile e$^{-x/\lambda }$, where $\lambda$ is a variable depth parameter.

The particles are propagated through the experimental geometry using the Livermore low energy physics models [51]. We use production cut lengths⁵ of 1 $\upmu$m for $e^{-}/e^{+}$ and 10 $\upmu$m for $\gamma$’s. For LMO these correspond to 1 keV energy thresholds for both $e^{-}/e^{+}$ and $\gamma$’s. This choice is based on a study of the impact of the production thresholds on the detected spectra. Thresholds of 1 keV and 250 eV for LMO give comparable spectra, while the computing time is significantly reduced for a 1 keV cut length.

5.1 Geometry

We implement a detailed geometry of the CUPID-Mo towers in the MC simulations. In particular, we reproduce the size of each LMO crystal [36] on an individual basis to take into account variations between the crystals. We also include the Ge LDs, the PTFE clamps, the reflective foils surrounding the crystals and the copper holders which are implemented as accurately as possible. Figure 2 shows the GEANT4 rendering of the simulation geometry of the 10 mK chamber, with the five towers of CUPID-Mo in the front. We included the readout cables, the springs, the EDELWEISS Ge detectors and their connectors. The copper structure supporting the crystals, composed of four copper plates and three copper columns made of NOSV copper, is also included in the simulated geometry. The four copper plates are held by brass screws with a relatively high mass (see Table 3) which have been modeled as well.

Figure 3 shows the simulated geometry of cryostat and electronics. The 10 mK, 1 K, 50 K, 100 K and 300 K thermal screens are included individually. The internal polyethylene shielding and lead shielding are also implemented in the geometry. We also note that the geometry extends and includes far components below the 1 K lead shield that are less radiopure, like the dilution unit, the 300 K electronics, the pumps and the He reservoir that is expected to be important for neutron simulations.

5.2 Detector response model

To compare the simulated spectra to the measured data we need to account for the finite energy resolution and response of the detectors. The following features are accounted for through a post processing of the MC simulation spectra:

• Energy resolution;

• Energy threshold of 40 keV;

• Event multiplicity;

• Scintillation light and LD resolution;

• Cut efficiencies;

• Inactive periods of detectors;

• Pile-up and delayed coincidences in decay chains.

We compute the energy resolution per detector-dataset pair as explained in Sect. 3.10. In particular, for a pulse with energy $E_{MC}$ in channel c and dataset d we sample from a Gaussian with mean $E_{MC}$ and standard deviation $R\times \sigma _{c,d}$. As is done in experimental data we discard pulses below the energy threshold, $<40$ keV, and compute the multiplicity as the number of detectors with $E>40$ keV for each simulated event.

We also reproduce the signals measured in the light detectors. We have parameterised the scintillation light energy measured by the LD in data as a function of LMO energy as a second order polynomial, for each LMO and side LD channel. We also parameterise the LD energy resolution as a Gaussian with standard deviation $\sqrt{p_0^2+p_1E}$. We use this parameterisation to generate a random scintillation light yield for each event which is summed with the energy deposited in MC from direct particle interactions. We use these light detector energies to reproduce the light yield cuts in the same way as in experimental data in Sect. 3.

To account for inactive detectors we assign a random timestamp from the data taking period to each simulated MC event. This allows us to apply the same cuts to the simulated data and remove events from detectors considered inactive and to account for the reduction of event multiplicity in these periods.

5.3 Simulated background sources

Some components produce indistinguishable spectra of energy deposits in the crystal, or in other words, they exhibit degenerate spectral shapes. In this case, we either group them, or generate the radioactive decays in only one, which accounts for all elements with degenerate spectra. This simplification reduces the number of free parameters in the fit of the simulations to the data, however, we need to keep in mind that the posterior distributions account for the sum of the grouped elements.

The reflective foils, the PTFE clamps, and all other passive elements directly facing the crystals produce degenerate spectra. For this reason we have generated radioactive decays only in the Reflectors. We simulated radioactive decays separately in the connectors, the cables and the springs in the detector chamber at 10 mK shown in Fig. 2, right. All the copper elements made of NOSV at the 10 mK stage (holders, four plates, three columns and 10 mK cryostat screen) have been grouped in one background source, which we refer to as Copper supports.

We have found that all thermal screens exhibit degenerate spectra, thus we group the screens made of NOSV copper and refer to as Cryostat screens. We have also found that the internal polyethylene shielding spectrum is degenerate with the one from internal lead shielding. Thus, we have chosen to include only the internal polyethylene shielding contribution in the fit. This element takes into account all contributions from background sources below the internal lead shielding, as the elements of the dilution unit or the 300 K electronics also shown in Fig. 3.

In addition we have considered as a source the outer cryostat screen, called screen 300 K. This volume also includes the contribution from radon present in the air between the 300 K screen and the external lead shielding.

6 Background model

The goal of the background model is to describe the data (Sect. 3) with the MC simulations (Sect. 5). The parameters of the model then tell us the radioactive contamination of various components of the experiment. We use a binned simultaneous maximum likelihood fit, performed in a Bayesian framework with a Markov Chain Monte Carlo (MCMC) approach [52], developed by CUORE and further optimized by CUPID-0 collaborations [15, 42] using the JAGS software [53, 54]. We model the data in spectra i (${\mathcal {M}}$ $_{1,\beta /\gamma }$, ${\mathcal {M}}$ $_2$ and ${\mathcal {M}}$ $_{1,\alpha })$ and energy bin b as:

$$\begin{aligned} f_i(E_b;{\vec {N}})=\sum _{j=1}^{N_s} N_{j}\times f_{j,i}(E_b). \end{aligned}$$

(4)

The sum j runs over the simulated MC sources, $N_{j}$ is a scaling factor for each source (shared by all three spectra) and $f_{j,i}(E_b)$ are the simulated MC spectral shapes, with $E_b$ the energy of bin b. The likelihood function for data $\mathcal {D}$, including the 3 spectra ${\mathcal {M}}$ $_{1,\beta /\gamma }$, ${\mathcal {M}}$ $_2$ and ${\mathcal {M}}$ $_{1,\alpha }$ is then given by the product of Poisson distributions, $\text {Poiss}(n_{i,b};f(E_i(E_b;{\vec {N}})))$, for $n_{i,b}$ observed counts in bin b of spectrum i, and prediction $f_i(E_b;{\vec {N}})$ for the set of parameters ${\vec {N}}$:

$$\begin{aligned}&\text {ln}{\mathcal {L}}\big ({\mathcal {D}}|f(E_b;{\vec {N}})\big )= \sum _{i=1}^{3} \sum _{b=1}^{N_b(i)} \ln {(\text {Poiss}(n_{i,b};f_i (E_b;{\vec {N}})))} \nonumber \\&\quad = \sum _{i=1}^{3} \sum _{b}^{N_b(i)} n_{i,b} \times \ln {(f_i(E_b;{\vec {N}}))} - f_i(E_b;{\vec {N}}) - \ln (n_{i,b}!). \end{aligned}$$

(5)

Here the sum i is over the three data spectra and b goes over the bins in each spectrum.

JAGS samples the full posterior probability distribution $p({\vec {N}}|{\mathcal {D}})$ given by Bayes theorem:

$$\begin{aligned} p({\vec {N}}|{\mathcal {D}}) = \frac{{\mathcal {L}}({\mathcal {D}}|{\vec {N}})\times \pi ({\vec {N}})}{p({\mathcal {D}})}, \end{aligned}$$

(6)

using MCMC. The prior probabilities, $\pi ({\vec {N}})$ are discussed in Sect. 6.3. For each parameter we also extract the marginalised posterior distribution by integrating over the parameter space $\Omega$ (excluding the parameter of interest):

$$\begin{aligned} p(N_j|{\mathcal {D}})=\int _{\Omega } p({\vec {N}}|{\mathcal {D}})d\vec {\Omega }. \end{aligned}$$

(7)

We choose the mode of the marginalised distribution as our point estimate of the parameter and we compute, by integrating, the smallest 68% Bayesian credible intervals, c.i., around the mode. If the lowest 68% includes the value zero, we give an upper limit at 90% c.i.

6.1 MC simulation of $^{56}$Co calibration source

We have performed a calibration with a $^{56}$Co source to validate the energy calibration and resolution of the CUPID-Mo detectors in the $0\nu \beta \beta$ ROI, at $\sim 3$ MeV. The measurement is also useful to test and validate the implementation of the Monte Carlo simulations. Two $^{56}$Co sources with an activity of $41 \pm 8$ Bq, measured with HPGe $\gamma$ spectroscopy immediately after the calibration, were placed on the outer cryostat screen, inside the external shielding. The configuration was chosen to achieve the highest counting rate in the ROI for all detectors with a total rate below 0.125 Hz as an upper limit on the tolerable pile-up.

We have performed a fit of the calibration data to a MC simulation of the $^{56}$Co sources summed with a background component (detailed later in this section) and pile-up, with only uniform priors. We describe in Sect. 6.3.1 how the spectral shape of the pile-up events is obtained. The fit has thus three parameters: the normalization factor of the background, the one of the $^{56}$Co sources, and the one of the pile-up events. We know the normalization of the background from the background model fit. Comparing it to the normalization factor of the background in the calibration data, we obtain the efficiency of the cuts in the calibration data $(68.7\pm 1.4)\%.$ From the normalization factor of the $^{56}$Co sources we obtain the activity of the sources without the efficiency correction. Using the estimated efficiency, we derive the final activity of the $^{56}$Co source of $(50 \pm 1)$ Bq, which is in good agreement with the measured activities.

The model shows good agreement with the data in the whole energy range of the fit 200–4000 keV, as shown in Fig. 6. In Fig. 7 we present the region above 2800 keV with 2 keV binning, where the comparison in this region can be better appreciated. This fit shows that the MC implementation of the set-up is accurate and that the MC is able to describe well the data.

Fig. 6

Top: Comparison between the $^{56}$Co calibration data and MC simulations, ${\mathcal {M}}$ $_1$ data, with a variable binning in the region between 200–4000 keV. Bottom: Bin by bin ratio of data and MC. Most of the values are within 3 $\sigma$, with discrepancies below 20%

Fig. 7

Comparison between the $^{56}$Co calibration data and MC simulations, ${\mathcal {M}}$ $_1$ data in the region of interest

6.2 Background sources list

The background model fit includes 41 background sources associated to the bulk volume of the components identified in Sect. 5.3. We included 2 additional sources of surface contamination: the LMO crystal and the Reflective foil (representing all sources facing the crystals). For a detailed list of the sources associated to radioactive contaminants we refer to Tables 4 and 5. The complete list of the background sources in our fit is:

• Crystal:

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

  • $2\nu \beta \beta$ decay of $^{100}$Mo to $^{100}$Ru 0$_1^+$ excited state,

  • pile-up (random coincidence of 2 events in the same crystal happening so close in time that the signal is equivalent to that of the sum of the two events),

  • $^{99}$Mo,

  • 12 bulk sources of natural radioactivity detailed in Table 4,

  • 8 sources associated to surface contamination, listed in Table 4,

  • $^{210}$Pb surface contamination with 1 nm and 1 $\mu$m implantation depth.

• Reflectors:

  • 3 sources associated to bulk contaminations, Table 5,

  • 8 sources associated to surface contaminations, Table 5,

  • $^{210}$Pb surface contamination with 100 nm and 1 $\upmu$m implantation depth.

• Close sources, 10 mK and infrastructure: 27 sources associated to the bulk volume of these components, listed in Table 5,

• Random coincidence of 2 events in two different crystals called accidentals.

Table 4 Radioactive contaminations of the LMO crystals derived from the background model of the CUPID-Mo data, with 2.71 kg $\times$ year exposure. The upper table shows the bulk activities. We report also the results under the assumption of no surface contaminations, to study the effect in the fit of the anticorrelation between bulk and surface activities. The lower table shows the surface activities, we give the activities with MC simulation with 10 nm implantation depth (see text). The effect of including a contribution with a depth parameter of 10 $\upmu$m is shown on the last column

A total of 67 sources are included in the fit. As mentioned in Sect. 5 we have modelled surface contaminations with an exponential density profile e$^{-x/\lambda }$. We have simulated surface contaminations with $\lambda = 10$ nm and 10 $\upmu$m for all radionuclides in the U and Th chains. The choice 10 nm is driven by the fact the typical range of recoiling nuclei is of the order of some nm for $\alpha$ decays in the U and Th chains. The choice of 10 $\upmu$m considers that mechanical crystal preparation including cutting, cleaning and polishing can lead to deep surface damage and implantation depths of $\sim \upmu$m, however depths $> 10~\upmu$m would be effectively equivalent to bulk contaminations. We observed that the component corresponding to the shallow surface contamination of the crystal, with $\lambda = 10$ nm, is needed to properly fit our data. Surface contaminations with $\lambda = 10~ \upmu$m give activities which are compatible with zero. Thus, for simplification to minimize the number of degrees of freedom, we choose to include in the fit only the crystal surface contaminations with $\lambda = 10$ nm implantation depth. Due to the small thickness (70 $\upmu$m) and low density of the Reflectors, surface contaminations with $\lambda = 10~\upmu$m are degenerate with bulk contaminations. Both produce continuous spectra due to the partial energy loss of $\alpha$ particles in the Reflector and the detection of the remaining kinetic energy in the crystal. We have therefore chosen to include only the surface contaminations with $\lambda =10$ nm for the Reflectors in the background model fit.

To simulate surface and bulk contaminations in the crystal and the Reflectors, we have generated the decay chains to take into account time correlations and exploit the delay coincidences (see Sect. 5), as done for the data. We observed that the bulk contamination in the Reflectors produce a flat spectra independent of the specific part of the radioactive decay chain and our fits showed the activities of the various subchains (excluding $^{210}$Pb to $^{206}$Pb) were compatible. We hence simplify the fit model by assuming that the entirety of the U/Th (excluding $^{210}$Pb to $^{206}$Pb) chain is in secular equilibrium for the Reflectors.

In addition to the U/Th chains, other contaminations have been included in the crystals. In particular, $^{40}$K can be found as a result of an initial contamination of the powder used to grow the crystals [31]. Some anthropogenic radionuclides due to fall out can also be found in enriched crystals [55], thus, we have included the bulk contaminations $^{87}$Rb and $^{90}$Sr+$^{90}$Y, which are pure $\beta$-emitters. $^{99}$Mo was produced by activation during a neutron calibration with an AmBe source.

For all sources that aren’t facing the crystals we have simulated decays of the daughters of $^{226}$Ra and $^{228}$Th. We identify in Table 3 a large content of $^{210}$Pb in the brass screws holding the detector plates. We have simulated $^{210}$Bi in this component and use it to account for this contamination in all 10 mK and infrastructure sources. Cobalt isotopes are expected to be primarily the result of cosmogenic activation in copper. We have therefore chosen to locate $^{60}$Co and $^{57}$Co in the Copper supports and use it to account for this contamination in all 10 mK and infrastructure sources. In all the components we use $^{228}$Ac $\gamma$ emitter with three main $\gamma$ peaks clearly visible in the data, not in equilibrium with $^{228}$Th. In doing so, we have observed that for the 300 K screen the two values from the fit are compatible with equilibrium thus, to reduce the number of components in the fit we have combined $^{228}$Ac and $^{228}$Th.

6.3 Choice of priors

We consider informative priors both from screening measurements (Sect. 4) and from other independent measurements. We have informative priors on the contribution of the:

• $2\nu \beta \beta$ $^{100}$Mo to $^{100}$Ru 0$_1^+$ excited state, which has been taken as $T_{1/2}=(6.7 \pm 0.5) \times 10^{20}$ year (average value from [9]),

• Stainless steel springs included in the set-up to mitigate vibrational noise. These are modelled with high accuracy in the MC, and the values measured by HPGe and used as priors are given in Table 2,

• Random coincidence (pile-up and accidentals) events, determined from the rate of single events and from a measurement with a calibration source (see below).

6.3.1 Random coincidence events

Energy deposition in either one or two crystals randomly coinciding in time can cause non-negligible contributions to the high energy region both in ${\mathcal {M}}$ $_{1,\beta /\gamma }$ and even more so in ${\mathcal {M}}_2$. This is a particular concern for $^{100}$Mo, due to the fast rate of $2\nu \beta \beta$ decay of $T_{1/2} \sim 7 \times 10^{18}$ year [56], or around 2 mHz per detector. The events in two different detectors are referred to as accidentals and contribute to ${\mathcal {M}}_2$ spectrum. The random coincidences in the same detector are referred to as pile-up and contribute to the ${\mathcal {M}}_1$ spectrum.

We predict the spectral shape of these events by convolving the experimental ${\mathcal {M}}$ $_{1,\beta /\gamma }$ spectrum with itself, i.e. by selecting randomly two energies in the experimental ${\mathcal {M}}$ $_{1,\beta /\gamma }$ spectrum and summing them. The ${\mathcal {M}}$ $_{1,\beta /\gamma }$ and the resulting random coincidences spectra are shown in Fig. 8.

The expected number of accidentals is then given by:

$$\begin{aligned} \hat{N}_{\text {acc}}=N^2 \frac{\Delta t}{t}\times \frac{N_{\text {LMO}}-1}{N_{\text {LMO}}}, \end{aligned}$$

(8)

where N is the total number of ${\mathcal {M}}_1$ events, $\Delta t/t$ is the ratio of the width of the coincidence time window, $\Delta t = \pm 10~\text {ms}$, to the total measurement time and $N_{\text {LMO}}$ is the number of LMO detectors. For the accidental random coincidences we place a prior as a Gaussian function with mean $\hat{N}_{\text {acc}}$ and $\sigma$ $\sqrt{\hat{N}_{\text {acc}}}$. We have used $N = 1.2 \times 10^6$, for a total measurement time of $2.2 \times 10^7$ s. We include this contribution only in the ${\mathcal {M}}_2$ spectrum.

The rate of pile-up events is generally lower than the rate of accidentals, but it is also less well constrained as the coincidence time, or effective time resolution, $\Delta t_{\text {eff}}$, is unknown a-priori and determined by the effectiveness of the pulse shape cuts used in the analysis. In general this time resolution will also be dependent on the energy of both the primary and secondary pulse as well as on their separation. However, since we are only interested in events in a narrow range of a high energy region $(\sim$ 3 MeV) we can treat this to a good approximation, as energy independent and simplify Eq. 8 to:

$$\begin{aligned} \hat{N}_{\text {pileup}}=N^2 \frac{\Delta t_{\text {eff}}}{t}. \end{aligned}$$

(9)

For thermal detectors typically the timing resolution is between the inverse sampling frequency and detector rise time. In CUPID-Mo the inverse of the sampling rate is 2 ms and the median value of the rise-time is 24 ms, with 8 ms spread [36]. Similar LMO crystals to CUPID-Mo, tested at LNGS have achieved 1–2 ms effective timing resolution [57] also using PSD only on the LMO channel (as we have done in CUPID-Mo). However, this test has pulses with a higher sampling frequency with respect to CUPID-Mo, different operation temperature and noise conditions and slightly different rise time, so cannot be directly extrapolated.

For the prior on pile-up events in our fit we use a measurement with Th/U calibration sources. We consider all events between 4–5.5 MeV, and from Eq. 9 we can obtain a value for $\Delta t_{\text {eff}}$. As there are no events in the selected region, we obtain a prior for $\Delta t_{\text {eff}}<7$ ms, at 90% c.i. As zero events are obtained, the corresponding probability density function is an exponential, that we use to place a prior on the pile-up rate. We include this contribution only in the ${\mathcal {M}}_{1,\beta /\gamma }$ spectrum.

Fig. 8

Experimental ${\mathcal {M}}$ $_{1,\beta /\gamma }$ and random coincidences (obtained by convolution of ${\mathcal {M}}$ $_{1,\beta /\gamma }$, arbitrary normalization) spectral shapes. We observe that the random coincidences distribution is shifted to higher energies (as expected [58, 59]) and could cause a background at the ROI

6.4 Binning and choice of energy intervals

The energy range of the fit is 100–4000 keV for ${\mathcal {M}}$ $_{1,\beta /\gamma }$, 400–4000 keV for ${\mathcal {M}}$ $_2$ and 3000–10000 keV for ${\mathcal {M}}$ $_{1,\alpha }$. We use a variable binning for the three spectra to have enough counts in each bin to minimize the effect of statistical fluctuations. We choose a minimum bin size of 15 keV for ${\mathcal {M}}$ $_{1,\beta /\gamma }$ and ${\mathcal {M}}$ $_2$, and 20 keV for ${\mathcal {M}}$ $_{1,\alpha }$. We set the minimum number of counts in each bin to be 50 for ${\mathcal {M}}$ $_{1,\beta /\gamma }$ and 30 for ${\mathcal {M}}$ $_2$ and ${\mathcal {M}}$ $_{1,\alpha }$. We choose each peak to be fully contained in one bin, to minimize the systematic effect of the detector response on our results.

7 Results

The result of the simultaneous fit of ${\mathcal {M}}$ $_{1,\beta /\gamma }$, ${\mathcal {M}}$ $_2$, and ${\mathcal {M}}$ $_{1,\alpha }$ to the CUPID-Mo data with 2.71 kg $\times$ year exposure is shown in Fig. 9. Tables 4 and 5 show the fit results, discussed in Sect. 7.2. We find that our background model is able to reconstruct well the 3 data spectra. On each spectrum the data over model ratio is shown, where the colors correspond respectively, to $\pm 1,$ $\pm 2,$ and $\pm 3$ $\sigma$ with:

$$\begin{aligned} \sigma _{i,b}= \frac{\sigma _{\text {data},\text {i},\text {b}}}{n_{\text {model},\text {i},\text {b}}}, \end{aligned}$$

(10)

where $n_{\text {model},\text {i},\text {b}}$ is the predicted number of counts in bin b and spectrum i and $\sigma _{\text {data},i,b}$ is the standard deviation of the data in this bin, $\sqrt{n_{\text {data},i,b}}$.

Table 5 Radioactive contaminations of the setup components derived from the posterior distribution of the background model fit. Uniform, non-informative priors are used except for the $^{228}$Th, $^{226}$Ra and $^{40}$K contaminants in the springs. For surface contaminations, the simulated depth is 10 nm. The last column shows the activities from screening measurements when available (see Tables 2 and 3 in Sect. 4)

Fig. 9

Experimental data and background model simultaneous fit reconstruction of the 3 CUPID-Mo data spectra. Two upper panels: ${\mathcal {M}}$ $_{1,\beta /\gamma }$, $\beta /\gamma$’s spectrum with energy deposits in only one detector. Middle panels: ${\mathcal {M}}$ $_2$, multiplicity 2 events, histogram of the 2 summed energies. Two bottom panels: ${\mathcal {M}}$ $_{1,\alpha }$, multiplicity 1 events in $\alpha$ energy region. For each one, the lower panel shows the ratio between experimental counts and reconstruction counts for each bin. The colors indicate the uncertainties at $\pm 1,$ $\pm 2,$ and $\pm 3$ $\sigma$

To investigate the goodness of the fit we generate pseudo-experiments, or toy Monte Carlo simulations. We sample randomly according to a Poisson distribution the events in each energy bin of the background model best fit reproduction. We fit independently each pseudo experiment and obtain the likelihood${\mathcal {L}}({\mathcal {D}}|{\vec {N}})$, i.e. the probability of the experimental data ${\mathcal {D}}$ given the set of parameters ${\vec {N}}$ of our model. We show in Fig. 10 the result for ${\mathcal {M}}$ $_{1,\beta /\gamma }$, ${\mathcal {M}}$ $_2$ and ${\mathcal {M}}$ $_{1,\alpha }$. The mean of the distributions of ${\mathcal {M}}$ $_{1,\beta /\gamma }$ and ${\mathcal {M}}$ $_2$ agree well with the value of the data. For ${\mathcal {M}}$ $_{1,\alpha }$ the result demonstrate a modest incompatibility between the data and the model probably arising from an incomplete modelling of $\alpha$ detector response or an $\alpha$ miscalibration. This effect is visible for E > 6 MeV in Fig. 9 bottom panels. This modest incompatibility has driven the choice of a systematic in our model and we have thus performed a fit with an energy range 3000–6360 keV for ${\mathcal {M}}$ $_{1,\alpha }$. We detail in Sect. 7.4 the results. The p value obtained are $p=0.38$, $p= 0.04$ and $p \sim 0$ for ${\mathcal {M}}$ $_{1,\beta /\gamma }$, ${\mathcal {M}}$ $_2$ and ${\mathcal {M}}$ $_{1,\alpha }$ respectively.

Fig. 10

Distribution of $- \ln {\mathcal {L}}({\mathcal {D}}|{\vec {N}})$ from the toys for the ${\mathcal {M}}$ $_{1,\beta /\gamma }$ (left), ${\mathcal {M}}$ $_2$ (middle) and ${\mathcal {M}}$ $_{1,\alpha }$ (right) spectra. The red line shows the $- \ln {\mathcal {L}}({\mathcal {D}}|{\vec {N}})$ of the reference fit for each of the spectra

7.1 SSD and HSD $2\nu \beta \beta$ decay mechanisms

The transition between the ground states of $^{100}$Mo and $^{100}$Ru, with spin parity 0$^{+}$, is realized via virtual $\beta$ transitions through 1$^{+}$ states of the intermediate nucleus $^{100}$Tc. Nuclear theory does not predict a-priori whether this transition is realized dominantly through the 1$^{+}$ ground state (SSD hypothesis) or through higher excited states of $^{100}$Tc (HSD hypothesis) [60].

We have found that the SSD mechanism of $2\nu \beta \beta$ decay to $^{100}$Ru ground state reproduces fairly well the data with a $p=0.38$, while the HSD model does not, $p \sim 0$. Since our data clearly favours SSD over HSD mechanism for $2\nu \beta \beta$, we have used the SSD model in our final fit.

7.2 Contaminations derived from the fit

7.2.1 LMO crystal contaminations

Fig. 11

Experimental ${\mathcal {M}}$ $_{1,\alpha }$ spectrum reconstruction showing the components of the ${\mathcal {M}}_{1,\alpha }$ background model fit. Crystal and Reflector contaminations include bulk and surface. The surface contaminations are modelled with an exponential density profile and $\lambda = 10$ nm parameter depth. The peaks in the spectrum are described by the radioimpurities in the crystal and the continuum by the ones in the bulk of the Reflectors. The small contribution from 10 mK sources corresponds to the holders

The ${\mathcal {M}}$ $_{1,\alpha }$ spectrum is populated by $\alpha$ decays occurring in the crystals and in the elements directly facing them. As described in Sect. 3.11 we included bulk and surface contaminations in the crystals in our fit. We show in Fig. 11, the resulting components. Since we do not observe clear $\alpha$-energy (NR escape) peaks due to the very low levels of contaminations and thus limited statistics, bulk and surface contaminations are anticorrelated. We performed studies concerning the effect of the location of the contaminations on bulk or surface in the fit results, which we discuss later.

The largest peak in the $\alpha$ region is the $^{210}$Po peak. This peak is largely described by the Q-value component of the crystal bulk. For the $^{210}$Po in order to fit as much as possible the particular shape peak in addition to 10 nm, implantation depths of 1 $\upmu$m and 1 nm are used in the crystals, and implantation depths of 100 nm and 1 $\upmu$m are used in the Reflectors. In Fig. 11 the left tail of the $^{210}$Po peak is described by the surface contamination on the Reflectors.

The summary of the crystal activities extracted from the fit is presented in Table 4. The LMO crystal contaminations by radionuclides from the $^{238}$U and $^{232}$Th chains are all below 1 $\upmu$Bq/kg. As a study of the effect of the bulk versus surface location, we performed a fit without surface contaminations. These results show that, even under this extreme assumption, the results on the bulk activities do not vary significantly. As shown in Fig. 5 the peak at 4.8 MeV contains $^{234}$U and $^{226}$Ra alpha decays. In this analysis this peak is ascribed to $^{234}$U, with a significant uncertainty (reported in Table 4) in the resulting contamination due to the anticorrelation with the $^{226}$Ra contribution. Additionally, in this peak we could have a contribution from the neutron capture in $^{6}$Li [31]. Neutrons captured in $^{6}$Li produce an alpha particle plus tritium, $^{6}$Li(n,$\alpha )$ $^{3}$H, with a total energy 4.782 MeV. We note also that the level of $^{228}$Ra is not constrained by any $\alpha$ peak.

There is clearly a larger $^{210}$Po contribution than the rest of the $^{238}$U chain, at the level of 96 $\upmu$Bq/kg, possibly introduced during the purification of the enriched material [61]. There are also traces of $^{190}$Pt, caused by the crystal growth in a platinum crucible [62] and we find $^{40}$K and $^{90}$Sr+$^{90}$Y at the level of some hundreds of $\upmu$Bq/kg. We note that $^{210}$Pb, $^{87}$Rb, $^{90}$Sr + $^{90}$Y and $^{40}$K do not represent a potential background for $0\nu \beta \beta$ search, as the $Q{_{\beta }}$ of these radioisotopes is much lower than the $0\nu \beta \beta$ ROI at 3 MeV.

We show in Table 4 (bottom) the surface contaminations of the crystals derived from the fit. We studied the effect of including also a contribution with a depth parameter of $10~\upmu$m (i.e., including surface contaminations with $\lambda = 10$ nm and 10 $\upmu$m) and the decay activity is shown in the third column of the table. The results are compatible with the fit including only 10 nm contributions. We observe clear anti-correlation for a given decay chain between the bulk and the surface contaminants in the crystal, but also with the surface of the Reflectors. These anti-correlations are taken into account in the uncertainties given in Table 4.

7.2.2 Radioactive contaminations of the setup components

A list of sources included in the fit and their resulting activities obtained from the marginalised mode and 68% c.i. are shown in Table 5.

The derived activities for the component called Reflectors are mainly constrained by the fit of the continuum in the $\alpha$ region. The values are larger than the measured radioactivities of the reflectors themselves, in particular, in $^{226}$Ra. We remind that this component takes into account all elements directly facing the crystals: PTFE, NTDs, LDs, bonding wires. A contamination of the reflecting foils introduced during the detector assembly could be conceived, explaining the activities obtained in the fit.

Concerning the surface activity on the reflecting foils, we performed a measurement with the BiPo-3 detector [41] which measures $^{214}$Bi and $^{208}$Tl levels through delay coincidences in the Bi-Po cascades. We can also convert the ICPMS results of the bulk measurement by assigning all the contamination to the surface. The surface activity of the Reflectors derived from the fit agrees well within uncertainties with both measurements.

The derived activities in the Kapton cables, the Connectors, the Brass Screws and the Copper supports agree well with the measured values. For the Cryostat Screens the activities obtained in the fit are higher than the measured levels from the raw copper. This points out to an additional contamination introduced during the fabrication of the screens for example due to the weldings. In particular, we have identified from the experimental data that the detectors facing the weldings in the cryostat screens have higher rates in the 2615 keV peak of $^{208}$Tl.

The Screen 300 K accounts for the residual environmental $\gamma$’s and the radon present in the gap between the outermost cryostat screen and the external lead shielding. The $^{226}$Ra contamination derived from the fit shown in Table 5 can be translated into a radon level concentration resulting in $(22\pm 3)$ mBq/m${^3}$, which is in good agreement with measurements of 20 mBq/m${^3}$ provided by the radon mitigation system in the LSM [34].

Figure 12 shows the breakdown of the components in the fit of ${\mathcal {M}}$ $_{1,\beta /\gamma }$. In the region 0.8–3 MeV the dominant contribution is the $2\nu \beta \beta$ from $^{100}$Mo and the most important contributions from the radioactivity in the materials are the cryostat and shields. We discuss in the next section the main sources in the $0\nu \beta \beta$ region.

Fig. 12

Background sources reconstructing the experimental ${\mathcal {M}}_{1,\beta /\gamma }$ spectrum, grouped by source location. In blue, $2\nu \beta \beta$ is the dominant contribution in [350–3000] keV. The most important contribution from the materials, below 3 MeV, are the cryostat and shields, shown in magenta

7.3 Background index in the $^{100}$Mo $0\nu \beta \beta$ ROI

We use our simultaneous fit to reconstruct the background index in the $0\nu \beta \beta$ region of interest. We chose to calculate the background index in the region $\pm 15$ keV around 3034 keV. This range is much wider than the experimental energy resolution, without including any $\gamma$ lines. We sample directly the full posterior distribution produced by JAGS for each step i in the Markov Chain by computing:

$$\begin{aligned} b_i =\sum _{\text {j}=1}^{{67}} \text {Pois}(N_j)\frac{w_{i,j}}{\Delta E \times Mt }. \end{aligned}$$

(11)

Here $b_i$ is the background index in the $0\nu \beta \beta$ region of interest, $N_j$ is the integral of the spectrum of MC source j in the ROI, $w_{i,j}$ is the weight for source j in step i, $\Delta E$ is the width of the ROI and Mt is the experimental exposure. The sum goes over all background sources. The MC simulations are themselves the result of a stochastic process they have a statistical uncertainty, this is accounted for by Poisson smearing the MC ROI integrals per step of the Markov Chain. We then use the distribution of $b_i$ for the full Markov Chain to estimate the marginalised posterior distribution of the background index. We perform this calculation for our maximal model with all parameters, we therefore naturally marginalise over all possible combinations of activities (for example surface or bulk radio-purity) consistent with our experimental data accounting for the systematic uncertainty due to source localisation.

From this calculation we extract the marginalised posterior of the background index shown in Fig. 13. This results in a measurement (mode ± smallest 68% interval) of:

$$\begin{aligned} b = 2.7^{+0.7}_{-0.6} \times 10^{-3}~\text {counts}/\text {keV}/\text {kg}/\text {year}. \end{aligned}$$

(12)

or, in terms of the number of moles of isotope, mol$_{\text {iso}}$, and energy resolution, $\Delta E_{\text {FWHM}}$:

$$\begin{aligned} {\mathcal {B}}= 3.7^{+0.9}_{-0.8}\times 10^{-3} \text {counts}/\Delta E_{\text {FWHM}}/\text {mol}_{\text {iso}}/\text {year}. \end{aligned}$$

(13)

This is the lowest background index achieved in a bolometric $0\nu \beta \beta$ decay experiment.

Next we reconstruct the contributions to the experimental background. We divide sources into five categories:

• Crystals U/Th;

• Pile-up;

• Reflectors;

• 10 mK sources;

• Cryostat and shields.

We emphasise that only the first three sources are relevant to CUPID. In the CUPID baseline the reflective foil is removed to improve background rejection. However, as it was noted before Reflectors include all the elements directly facing the crystals, PTFE, bonding wires, heaters. These elements will remain in CUPID. The final two are caused by materials in the EDELWEISS cryostat which is optimised for a dark matter rather than $0\nu \beta \beta$ decay search. The posterior distributions of background index from each source are shown in Fig. 14. We derive the background index for each of the sources in the same way as for the full posterior. We find that the crystals give the smallest contribution, with a background index:

$$\begin{aligned} 8.1^{+3.5}_{-2.5} \times 10^{-5}~\text {counts}/\text {keV}/\text {kg}/\text {year}. \end{aligned}$$

(14)

As shown in Fig. 14, the posterior probability for pile-up allows us to set an upper limit for its background index, $<1.4 \times 10^{-3}$ counts/keV/kg/year (90% c.i.). This is potentially the main background contribution, in particular due to the low CUPID-Mo sampling frequency (500 Hz) and lack of optimised cuts to remove pile-up. In CUPID, heat and light signals will be exploited together with optimised algorithms to remove pile-up events (see for example [63]). Figure 15 gives the background index extracted from Fig. 14 for each of the grouped components. They are obtained from the mode and the smallest 68.3% interval. For the pile-up the smallest 68.3% interval is compatible with zero, thus an upper limit is presented.

7.4 Systematics

To check the stability of the model and the systematic uncertainties, we perform a series of different fits. To take into account the systematic uncertainty due to MC statistics, we add a nuisance parameter in Eq. 5:

$$\begin{aligned}&\text {ln}({\mathcal {L}}(D|({\vec {N}}))) \nonumber \\&\quad =\sum _{i=1}^{3} \sum _{b=1}^{N_b(i)} \ln {(\text {Poiss}(n_{i,b};f_i(E_b;{\vec {N}})))} \nonumber \\&\qquad +\ln {(\text {Poiss}(N^{MC}_{j,i,b}; \hat{N^{MC}_{j,i,b}}))}, \end{aligned}$$

(15)

where, $N^{MC}_{j,i,b}$ is the number of MC events in bin b of source j in spectra i, and $\hat{N^{MC}_{i,b}}$ is the expected number. These nuisance parameters added to the model account for the integer Poisson fluctuations in the MC. We find that the fit remains largely unchanged with only a small change in the value of the background index.

Fig. 13

Posterior distribution of background index, showing the mode and the smallest 68.3% c.i., $2.7_{-0.6}^{+0.7}\times 10^{-3}$ counts/keV/kg/year

Fig. 14

Posterior distributions of background index of the several background sources grouped by source location. Also shown is the full posterior distribution

To check the stability of the fit, we perform different fits varying the binning, the energy fit region, the choice of background sources and, in particular, the bulk and surface contaminations in the crystals, as follow:

• Binning: we repeat the fit with 1, 2 and 20 keV fixed binning in ${{\mathcal {M}}}_{1,\beta /\gamma }$ and ${\mathcal {M}}$ $_2$. In all cases, the overall goodness of the fit remains, and the value of the background index is compatible within uncertainties to that of the reference fit, as shown in Table 6. We did not repeat the fit with 1, 2 and 20 keV on ${\mathcal {M}}$ $_{1,\alpha }$ due to the low statistics in each bin of the data;

• Fit energy region: our reference fit extends from 100 keV to 4 MeV for ${\mathcal {M}}$ $_{1,\beta /\gamma }$ spectrum. We vary the energy threshold to 200 keV and find that the background index only varies slightly;

• Choice of background sources: our calculation of the background index is naturally marginalising over this uncertainty (see Sect. 7.3). However, as an additional check we perform the fit without including the crystal bulk contribution for the U and Th chains. The values of the activities of the sources change, but the goodness of the fit remains very similar and the value of the background index remains almost unchanged. We then remove the crystal surface contamination and still obtain a background index compatible within uncertainties to that of the reference fit;

• Energy region of ${\mathcal {M}}$ $_{1,\alpha }$ fit: our reference fit extends from [3000–10000] keV. As described at the beginning of Sect. 7 the ${\mathcal {M}}$ $_{1,\alpha }$ fit shows a modest incompatibility between the data and the model, mainly in the region E > 6 MeV. We thus performed a fit in [3000–6360] keV to account for this incompatibility as a systematic uncertainty in our model. In doing so, the U and Th contributions in the crystal get more degenerated, resulting in an increase of the Th contamination assigned in the fit. Still the background index is compatible, within uncertainties, to that of the reference fit.

The results of these tests are summarized in Table 6. As argued above, the result given in Eq. 12 is naturally marginalising over the uncertainty on the choice of the background sources. Considering all tests in Table 6 as a systematic uncertainty (with 2 keV fixed binning) and adding them in quadrature, the background index in (3034 ± 15) keV results in:

$$\begin{aligned} b = 2.7^{+0.7}_{-0.6}(\text {stat})^{+1.1}_{-0.5}(\text {syst}) \times 10^{-3}~\text {counts}/\text {keV}/\text {kg}/\text {year}, \end{aligned}$$

(16)

or:

$$\begin{aligned} {\mathcal {B}}{} & {} = 3.7^{+0.9}_{-0.8}(\text {stat})^{+1.5}_{-0.7}\nonumber \\{} & {} \quad ~(\text {syst})~\times 10^{-3} \text {counts}/\Delta E_{\text {FWHM}}/\text {mol}_{\text {iso}}/\text {year}. \end{aligned}$$

(17)

Fig. 15

Background index for the various groups of sources. The values are extracted from the mode of each distribution of Fig. 14, with their respective uncertainties. The green bars correspond to the smallest 68.3% interval around the mode, and the yellow bars to the smallest 90% interval around the mode. For the pile-up the distribution is compatible with zero, thus we give an upper limit to 68.3% c.i. in green and 90% c.i. in yellow

Table 6 Background Index in ROI for different fits. The tests allow to check the stability of the model and assess the systematic uncertainties

We also verified that the reconstructed background index is not biased, by comparing the distribution of background indexes in toy Monte Carlo simulations to that of the reference fit.

7.5 Residual alpha background

Due to our $\alpha$ particle rejection, background events in the ROI from $^{226}$Ra and $^{228}$Th subchains in the bulk and the surface of the crystals generally arise only from energy depositions of $\beta$ or $\gamma$ particles. However, $^{238}$U, $^{234}$U, $^{230}$Th, $^{210}$Po and $^{232}$Th could also produce events in the ROI through energy deposits of $\alpha$ particles. Even if we apply a light yield cut to remove $\alpha$ background, it is still possible that some $\alpha$ events pass this selection cut.

From the background index distribution of the crystals, one can separate the background from $\beta /\gamma$ decays from that coming from $\alpha$ decays, as shown in Fig. 16. One can observe that a non-negligible part of the crystal background index is coming from $\alpha$’s that passes the light yield cut. This $\alpha$ background is coming from surface contamination of the crystals. It corresponds to an $\alpha$ particle that deposits energy in the crystal and where the nuclear recoil deposits energy in the LD. This kind of events can pass the light yield cut mainly for the crystals that face only one LD. We show in Fig. 17, left, the experimental ${\mathcal {M}}$ $_{1,\beta /\gamma }$ spectrum including all crystals, and the resulting spectrum selecting only the crystals that face two LDs. Such cut remove all the remaining alphas around 5.8 MeV. This effect is also visible in the background model. Figure 17, right, shows the reconstruction of the crystal component of ${\mathcal {M}}$ $_{1,\beta /\gamma }$ spectrum, and the resulting spectrum selecting only the crystals that face two LDs. We remind that in CUPID-Mo 5 of the 20 LMOs are facing only one LD due to being in the bottom floor of the towers, one LD was not operational and one had a poor performance affecting a further 4 LMOs. We expect that in the case where all the crystals face two LDs, as in CUPID, the $\alpha$ background contribution should be negligible.

Fig. 16

Posterior distribution of background index of the crystal from $\alpha$ and $\beta /\gamma$ contamination

Fig. 17

Left: Experimental ${\mathcal {M}}$ $_{1,\beta /\gamma }$ spectrum (in blue) adding a cut to select crystals that face two LDs (in red). Right: Fit reconstruction of the crystal component from ${\mathcal {M}}$ $_{1,\beta /\gamma }$ spectrum (in blue), adding a cut to select crystals that face two LDs (in red)

8 Conclusion

In this work we present the development of a background model capable of describing very accurately the CUPID-Mo experimental data with 2.71 kg $\times$ year exposure. We have performed a simultaneous fit of three data spectra, ${\mathcal {M}}$ $_{1,\beta /\gamma }$, ${\mathcal {M}}$ $_2$ and ${\mathcal {M}}$ $_{1,\alpha }$, to detailed Monte Carlo simulations. The model is performed in a Bayesian framework with a MCMC approach. We have shown by a fit to a $^{56}$Co calibration source that the MC implementation is accurate and that the MC is able to describe well the data. We used a total of 67 background sources including the bulk and surface radioactive contaminations in the crystal and the components of the set-up. We have performed systematic checks varying the binning, the energy fit region and the choice of background sources that showed the stability of the model.

We have found that the radiopurity of the Li$_{2}$ $^{100}$MoO$_4$ crystals is sufficient to reach the goals of the future $0\nu \beta \beta$ experiment CUPID. The radiopurity levels of $^{226}$Ra and $^{228}$Th are below 0.5 $\mu$Bq/kg. We obtain a background index in the region of interest of 3.7$^{+0.9}_{-0.8}$ (stat)$^{+1.5}_{-0.7}$ (syst) $\times ~10 ^{-3}$ counts/$\Delta E_{\text {FWHM}}$/mol$_{\text {iso}}$/year, the lowest in a bolometric $0\nu \beta \beta$ decay experiment.

The detailing of the background achieved in this work enables promising further studies. We can obtain the $2\nu \beta \beta$ decay rate of $^{100}$Mo with high precision. It also allows for studies on various process which could distort the spectral shape, like Bosonic neutrinos, CP violation or $0\nu \beta \beta$ with Majoron(s) emission.