The background model of the CUPID-Mo 0νββ experiment

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ The background model of the CUPID-Mo 0νββ experiment © 2023 the Authors Published version Augier, C.; Barabash, A. S.; Bellini, F.; Benato, G.; Beretta, M.; Bergé, L.; Billard, J.; Borovlev, Yu. A.; Cardani, L.; Casali, N.; Cazes, A.; Celi, E.; Chapellier, M.; Chiesa, D.; Dafinei, I.; Danevich, F. A.; De Jesus, M.; de Marcillac, P.; Dixon, T.; Dumoulin, L.; Eitel, K.; Ferri, F.; Fujikawa, B. K.; Gascon, J.; Gironi, L.; Giuliani, A.; Grigorieva, V. D.; Gros, M.; Helis, D. L.; Huang, H. Z.; Huang, R.; Imbert, L.; Johnston, J.; Juillard, A.; Khalife, H.; Kleifges, M.; Kobychev V., V.; Kolomensky, Yu. G.; Konovalov, S. I.; Kotila, J.; Loaiza, P.; Ma, L.; Makarov, E. P.; Mariam, R.; Marini, L.; Marnieros, S.; Navick, X.-F.; Nones, C.; Norman, E. B.; Olivieri, E.; Ouellet, J. L.; Pagnanini, L.; Pattavina, L.; Paul, B.; Pavan, M.; Peng, H.; Pessina, G.; Pirro, S.; Poda, D. V.; Polischuk, O. G.; Pozzi, S.; Previtali, E.; Redon, Th.; Rojas, A.; Rozov, S.; Sanglard, V.; Scarpaci, J. A.; Schmidt, B.; Shen, Y.; Shlegel, V. N.; Singh, V.; Tomei, C.; Tretyak, V. I.; Umatov, V. I.; Vagneron, L.; Velázquez, M.; Welliver, B.; Winslow, L.; Xue, M.; Yakushev, E.; Zarytskyy, M.; Zolotarova A., S. Augier, C., Barabash, A. S., Bellini, F., Benato, G., Beretta, M., Bergé, L., Billard, J., Borovlev, Y. A., Cardani, L., Casali, N., Cazes, A., Celi, E., Chapellier, M., Chiesa, D., Dafinei, I., Danevich, F. A., De Jesus, M., de Marcillac, P., Dixon, T., . . . Zolotarova A., S. (2023). The background model of the CUPID-Mo 0νββ experiment. European Physical Journal C, 83, Article 675. https://doi.org/10.1140/epjc/s10052-023-11830-2 2023 Eur. Phys. J. C (2023) 83:675 https://doi.org/10.1140/epjc/s10052-023-11830-2 Regular Article - Experimental Physics The background model of the CUPID-Mo 0νββ experiment C. Augier1, A. S. Barabash2, F. Bellini3,4, G. Benato5,6, M. Beretta7, L. Bergé8,J.Billard 1, Yu. A. Borovlev9, L. Cardani4,N.Casali 4, A. Cazes1,E.Celi 5,6, M. Chapellier8, D. Chiesa10,11, I. Dafinei4, F. A. Danevich4,12, M. De Jesus1, P. de Marcillac8, T. Dixon8, L. Dumoulin8,K.Eitel 13, F. Ferri15, B. K. Fujikawa14, J. Gascon1, L. Gironi10,11, A. Giuliani8, V. D. Grigorieva9,M.Gros 15,D.L.Helis 5,6, H. Z. Huang16, R. Huang7, L. Imbert8,a, J. Johnston17,A.Juillard 1, H. Khalife15, M. Kleifges18, V. V. Kobychev12, Yu. G. Kolomensky7,14,S.I.Konovalov 2, J. Kotila19,20,21, P. Loaiza8,L.Ma 16, E. P. Makarov9, R. Mariam8, L. Marini5,7, S. Marnieros8, X.-F. Navick15, C. Nones15, E.B. Norman22, E. Olivieri8,J.L.Ouellet 17, L. Pagnanini5,6, L. Pattavina6,23,B.Paul 15,M.Pavan 10,11, H. Peng24, G. Pessina11,S.Pirro 6,D.V.Poda 8, O. G. Polischuk4,12, S. Pozzi11,E.Previtali 10,11, Th. Redon8, A. Rojas25, S. Rozov26, V. Sanglard1, J. A. Scarpaci8, B. Schmidt15, Y. Shen16, V. N. Shlegel9, V. Singh7, C. Tomei4, V. I. Tretyak6,12,V.I.Umatov 2, L. Vagneron1, M. Velázquez27, B. Welliver7,L.Winslow 17,M.Xue 24, E. Yakushev26, M. Zarytskyy12, A. S. Zolotarova15 1Université Lyon 1, CNRS/IN2P3, IP2I-Lyon, 69622 Villeurbanne, France 2National Research Centre Kurchatov Institute, Kurchatov Complex of Theoretical and Experimental Physics, 117218 Moscow, Russia 3Dipartimento di Fisica, Sapienza Università di Roma, P.le Aldo Moro 2, 00185 Rome, Italy 4INFN, Sezione di Roma Tor Vergata, Via della Ricerca Scientifica 1, 00133 Rome, Italy 5Gran Sasso Science Institute, 67100 L’Aquila, Italy 6INFN, Laboratori Nazionali del Gran Sasso, 67100 Assergi, AQ, Italy 7Department of Physics, University of California, Berkeley, CA 94720, USA 8Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France 9Nikolaev Institute of Inorganic Chemistry, 630090 Novosibirsk, Russia 10 Dipartimento di Fisica, Università di Milano-Bicocca, 20126 Milan, Italy 11 INFN, Sezione di Milano-Bicocca, 20126 Milan, Italy 12 Institute for Nuclear Research of NASU, Kyiv 03028, Ukraine 13 Institute for Astroparticle Physics, Karlsruhe Institute of Technology, 76021 Karlsruhe, Germany 14 Nuclear Science Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA 15 IRFU, CEA, Université Paris-Saclay, 91191 Gif-sur-Yvette, France 16 Key Laboratory of Nuclear Physics and Ion-beam Application (MOE), Fudan University, Shanghai 200433, People’s Republic of China 17 Massachusetts Institute of Technology, Cambridge, MA 02139, USA 18 Institute for Data Processing and Electronics, Karlsruhe Institute of Technology, 76021 Karlsruhe, Germany 19 Department of Physics, University of Jyväskylä, PO Box 35, 40014 Jyvaskyla, Finland 20 Finnish Institute for Educational Research, University of Jyväskylä, P.O. Box 35, 40014 Jyvaskyla, Finland 21 Center for Theoretical Physics, Sloane Physics Laboratory, Yale University, New Haven, CT 06520-8120, USA 22 Department of Nuclear Engineering, University of California, Berkeley, CA 94720, USA 23 Physik Department, Technische Universität München, 85748 Garching, Germany 24 Department of Modern Physics, University of Science and Technology of China, Hefei 230027, People’s Republic of China 25 LSM, Laboratoire Souterrain de Modane, 73500 Modane, France 26 Laboratory of Nuclear Problems, JINR, 141980 Dubna, Moscow Region, Russia 27 Université Grenoble Alpes, CNRS, Grenoble INP, SIMAP, 38402 Saint Martin d’Héres, France Received: 2 May 2023 / Accepted: 10 July 2023 / Published online: 28 July 2023 © The Author(s) 2023 Abstract CUPID-Mo, located in the Laboratoire Souterrain de Modane (France), was a demonstrator for the next generation 0νββ decay experiment, CUPID. It consisted of an array of 20 enriched Li2100MoO4bolometers and 20 Ge light detectors and has demonstrated that the technolae-mail: [email protected] (corresponding author) ogy 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νββ decay and other processes with high precision. We also mea123 675 Page 2 of 24 Eur. Phys. J. C (2023) 83 :675 sure the radio-purity of the Li2100MoO4crystals 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) ×10−3counts/EFWHM/moliso/year, the lowest in a bolometric 0νββ decay experiment. Contents 1 Introduction ..................... 2 2 The CUPID-Mo experiment ............. 3 3 Experimental data .................. 3 4 Background sources ................. 8 5 Monte Carlo simulations ............... 8 6 Background model .................. 11 7 Results ........................ 14 8 Conclusion ...................... 22 References ........................ 24 1 Introduction Neutrinoless double beta decay (0νββ) 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νββ signal is a peak in the summed electron energy spectrum positioned at the Qββ (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νββ half-life are of the order of 1024–1026 year [2–8]. The Standard Model process, two-neutrino double beta decay, 2νββ, includes also the emission of two ¯νeand conserves lepton number. Unlike 0νββ decay, 2νββ has a continuous energy spectrum and has been observed in more than ten nuclei with half-lives in the range of 1018–1024 year [9]. One of the largest challenges in 0νββ 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νββ signal reducing the experimental sensitivity. During the last 10 years the scintillating bolometer technology has proved that bolometers based on lithium molybdate (Li2MoO4), are very promising detectors for next generation 0νββ searches [10,11]. Scintillating bolometers were developed to reduce the background observed in the current leading 0νββ bolometric experiment, CUORE [12]. In CUORE, the background in the region of interest is dominated by surface α’s emitted from the copper structure holding the detectors [13]. The array of 988 TeO2bolometers, installed at the Laboratori Nazionali del Gran Sasso, Italy, has observed a background in the 130Te region of interest (Qββ = 2527 keV) of (1.49 ±0.04)× 10−2counts/keV/kg/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 β/γ from αparticles with scintillating bolometers relies on the fact that the light emitted in the Li2100MoO4by αparticles is about a factor 5 smaller compared to the light emitted by β/γ’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 Qββ energy value of 100Mo (3034 keV) is higher than 2615 keV, implying a signal located above the majority of γ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 Li2100MoO4 (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 αbackground. LMO radiopurity for the U and Th chains of less than 10 and 3 µ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 α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νββ decay of 100Mo 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νββ decay rate and allows to disentangle between the Single State Dominance (SSD) or High State Dominance (HSD) mechanisms [19,20]ofthe2νββ decay process in 100Mo. It also provides the basis to study new physics processes outside the Standard Model, which could distort the spectral shape of the 2νββ spectrum, such as 0νββ decay with emission of Majoron(s), 2νββ decay with emission of Bosonic neutrinos, Lorentz invariance violation or sterile neutrinos [21–30]. 123 Eur. Phys. J. C (2023) 83 :675 Page 3 of 24 675 Fig. 1 Left: An individual CUPID-Mo bolometer showing the transparent Li2100MoO4crystal, 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] 2 The CUPID-Mo experiment CUPID-Mo was installed in the EDELWEISS cryogenic setup [31] at LSM. The experiment was in operation between March 2019 and June 2020. CUPID-Mo used 100Mo-enriched LMO crystals, where the 100Mo, the double beta isotope, has been enriched at ∼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 (3M Vikuiti®) is installed around the crystals, inside the copper holder, to increase the light collection. The average weight of the CUPID-Mo crystals is 210g and the total mass is 4.158 kg corresponding to 2.264 kg of 100Mo. 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 CUPIDMo towers in Fig. 2, are equipped with Kapton®pads and MillMax®connectors. The detector chamber consists of four copper plates made of NOSV®grade copper1to 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 1Copper of 99.9975% purity, produced by Aurubis, Hamburg, Germany. of NOSV and CUC22grade copper. An internal polyethylene (PE) shield, to shield against neutrons produced in the set-up components by (α, n) reactions or induced by muons [32], is mounted between the detectors and the internal lead shield, and has a temperature of ∼1 K. An internal lead shield of 14cm 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 20cm thick lead, with the innermost 2cm 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/m3[34]. Following the external lead shield, a 50cm 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±0.4)keV FWHM at 3034 keV [6] and α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]. 3 Experimental data The aim of our data processing is to convert the raw data stream into three calibrated energy spectra: β/γ like events with energy deposits in a single crystal (M1,β/γ ), events with energy deposits in two crystals (M2)and of α-like events (M1,α). 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 ×year of LMO corresponding to 2Copper of >99.990 % purity. 123 675 Page 4 of 24 Eur. Phys. J. C (2023) 83 :675 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 Fig. 3 GEANT4 rendering of the CUPID-Mo Monte Carlo simulation geometry, showing the cryogenic set-up 1.47 kg ×year of 100Mo. Our data is acquired as a continuous time-stream and digitized at 500 Hz by the EDELFig. 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 WEISS DAQ [36] and stored at both CC-IN2P3 (France) and NERSC (USA) for offline analysis. We acquire runs, periods of around 10–100h 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 232Th/238U source to calibrate the LMO detectors and a high activity 60Co source, which generates 100Mo X-rays in the detectors, to calibrate the LDs. We divide the data into twelve periods of ∼1 month of stable data taking, called datasets. 123 Eur. Phys. J. C (2023) 83 :675 Page 5 of 24 675 We discard three short periods of data (≈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 3s 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 232Th/238U calibration measurements. This calibration is accurate to around <1keV[6] which is sufficient for the binned background model fits. The LDs are calibrated using the dedicated 60Co calibrations which produces ∼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νββ decay or 2νββ decay to the ground state are very likely (> 75% probability [11]) to deposit energy in only one crystal. However, background events in particular from γ’s are likely to deposit energy in multiple crystals simultaneously, for example due to Compton scattering in one crystal, or multiple γ’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 ±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 3s 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 variables3which 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 αfrom β/γ particles. The use of light detectors also allows us to remove background events in which a particle deposits energy on our LDs. We select β/γ 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, ELD, and the mean expected β/γ LD energy L, normalized by the light band width (σ). We compute nfor M1events. As each LD has different characteristics, the calculation is done for each channel (crystal) c, each dataset dand for both side channels s, i.e.: nc,s,d=ELD −Lc,d,s(E) σc,d,s(E),(1) with Ethe measured LMO energy. The parameter nc,s,dhas a distribution expected to be centered at zero for β/γ’s and at a value different from zero for α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 nc,1against nc,2to be a bivariate Gaussian. As we observe no clear correlation we place a radial cut on the variable: D=n2 c,1+n2 c,2.(2) If only one LD is available the cut is instead placed just on this nc,s,d. We chose a cut of D<4 to select β/γ events and call this data spectrum M1,β/γ . We also construct a spectrum of M1,α events comprised of high energy M1events, E> 3 MeV, with no light selection applied. This data comprises almost entirely αparticles. The same events are obtained with a selection cut D>4, thus, for simplification we have chosen only the energy cut to select αevents. Unlike most other analysis of scintillating bolometers we also develop a light selection cut for M2events as described in detail in [18]. For a M2event 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 3The optimal filter test values or the χ2for rising and falling edges and for the pulse baseline. 123 675 Page 6 of 24 Eur. Phys. J. C (2023) 83 :675 pulse in a M>1 event. We require that each normalised LD energy (for each channel and side channel) is between −10 and 10 σ, for all but one LD. In this particular LD we observe an accidental contamination of 60Co. Therefore we generally observe γevents in the LMO and β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 to3σ, 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 60Co source we make a global LD anti-coincidence cut to remove the γbackground originating from this LD. We remove any events (on non-adjacent LMOs) with a trigger on this LD with energy >2keVwithina5mswindow. 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 M1,β/γ ,M2and M1,α 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 O(90%) tagging efficiency of muons with a minimal impact on the β/γ 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 232Th and 238U 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 214Bi (238U chain) and 208Tl (232Th chain). For 208Tl we veto events in 10×T1/2(1830s) following a suspected 212Bi α-decay. This time window contains >99.9% of the 208Tl decays. The very low CUPID-Mo radioactivity also enables a novel delayed coincidence cut removing 214Bi candidate events. The 222Rn decay chain proceeds as follows: 222Rn 3.8day ⇒ α5590 keV 218 Po 3.1min ⇒ α6115 keV 214 Pb 27.1min ⇒ β−1018 keV 214Bi 19.7min ⇒ β−3269 keV 214Po.(3) We can therefore tag the event based on the 222Rn or 218Po αevents and a fairly long dead time. We use energy cuts of 5985–6145 keV for 218Po and 5460–5620 keV for 222Rn to tag αcandidates. For either type of α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 13860s for 222Rn and 13620 s for 218Po. The two possible cuts on 222Rn or 218Po 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: –M1,β/γ : Events in one detector identified as β/γ, –M2: Events in coincidence between 2 crystals, the two energies deposited in each crystal are summed, –M1,α: Events in one detector with alpha energy scale (> 3MeV). Because of the relatively fast half-life of 2νββ in 100Mo (∼7×1018 year) and extremely low levels of contamination, relatively few peaks are observed in the M1,β/γ spectrum, where the spectrum of 2νββ decays of 100Mo is the dominant feature. The secondary datasets, M2and M1,α however contain a lower fraction of 2νββ events and therefore provide useful information to determine the location of radioactive contaminations. The experimental spectra after all cuts are shown in Fig. 5. 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 γpeaks from M2and M1,β/γ 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 γ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 210Po peak. We evaluate the pile-up efficiency, the probability a pulse will be superimposed with another in a 3s 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. 3.10 Energy scale and resolution We use the observed γ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. 123 Eur. Phys. J. C (2023) 83 :675 Page 7 of 24 675 Fig. 5 CUPID-Mo experimental data. Left: M1,β/γ : events in one detector identified as γ/β.M2: events in coincidence between 2 crystals, the two energies deposited in the crystals are summed. Right: M1,α : events in one detector with αenergy scale Table 1 Efficiencies for the cuts used on CUPID-Mo data. The PSD and Light Distance cut efficiencies are evaluated using γpeaks, and show no energy dependence in the range of the fit Cut Evaluation method Efficiency [%] PSD (M1,β/γ )M1,β/γ γ-peaks 95.2±0.5 PSD (M2)M2γ-peaks 96.9±0.5 Light distance (M1,β/γ )M1,β/γ γ-peaks 99.4±0.4 Light distance (M2)M2γ-peaks 97.7±1.8 Multiplicity 210Po 99.55 ±0.07 Rise time cut M1,β/γ γ-peaks 99.8±0.2 LD anti-coincidence 210Po 99.976 ±0.017 Muon veto cut 210Po 99.62 ±0.07 Delayed coincidences 210Po 99.16 ±0.01 Pile-up Noise 95.7±1.0 Total M1,β/γ 88.9±1.1 Total M283.3±2.5 Total M1,α 94.7±1.0 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 Eand width R×σc,d, where σc,dis the energy resolution in channel cand dataset d. This energy calibration is discussed in detail in [6]. 3.11 Features of data spectra We observe in Fig. 5that the spectrum of 2νββ decays of 100Mo dominates the M1,β/γ data, whereas the M2spectrum has significant contributions from natural radioactivity, shown by prominent γpeaks. These consist predominantly of decays from the 238U and 232Th decay chains, however we also observe contributions from 40K and cosmogenic activation products 60Co and 57Co. We also observe a short lived peak of 99Mo, present for ∼1 dataset, from neutron activation after a calibration with an AmBe neutron source. The spectrum M1,α is dominated by αdecays from components very close to the detectors. As shown in Fig. 5,in our data we observe a large contribution of 210Po, Eα= 5303 keV, with both a large Q-value and α-energy peak and peaks from several other nuclides in the U/Th chains. During αdecay the energy released is shared between the α-particle and recoiling nucleus (NR), with energy O(100 keV). In LMO crystals the range of αparticles is about 10 µ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 O(nm) in the crystal the αor recoil could escape, or both could be contained in the crystal. We therefore expect peaks at the NR energy, at the αenergy, and at the Q-value, with a relatively low flat continuum from partial contained α’s or NR. The ratio of these peaks depends on the depth of radioactive contamination. For a deeper contribution O(µm) the NR is almost always contained but the α’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 α-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. 5that we generally do not observe clear α-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 α-energy peaks creates a challenge for assessing the surface contamination. 123 675 Page 8 of 24 Eur. Phys. J. C (2023) 83 :675 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 238U and 232Th decay chains, 40K, and cosmogenic radionuclides have been measured by High Purity Ge γ-ray spectroscopy and ICPMS (Inductively Coupled Plasma Mass Spectrometry). The CUPID-Mo materials were chosen to minimize the 226Ra and 228Th contaminations, as the most critical radioactive backgrounds in the 3 MeV region relevant to 0νββ decay searches arise from 214Bi and 208Tl decays. Table 2reports 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, O(100 keV), which is the region of interest in dark matter searches. Table 3 shows the radioactivity in the EDELWEISS cryostat materials.4The 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. 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. 4In Table 3, the Kapton connectors, MillMax connectors and Cu Kapton cables belong to the EDELWEISS readout system, while the NOMEX cables are used for the CUPID-Mo readout. 4.1 Other contributions – muons, neutrons and environmental gammas The muon flux at the LSM is 5 muons/m2/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νββ 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±0.05 (stat.)±0.27 (syst.))× 10−6neutrons/s/cm2[43] and the ambient neutron flux (fast plus thermal) has been estimated ∼10−5neutrons/s/cm2in [43,44]. Previous work [45] showed that 48 cm of polyethylene reduces the neutron flux by a factor 2 ×106. 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±0.2(stat.)×10−2 γ/s/cm2[46]. Considering that 20cm of lead reduce the flux by about a factor 104, 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νββ 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 226Ra to 210Pb. 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 123 Eur. Phys. J. C (2023) 83 :675 Page 15 of 24 675 Fig. 8 Experimental M1,β/γ and random coincidences (obtained by convolution of M1,β/γ , 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 shown in Fig. 9. Tables 4and 5show 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 ±1,±2,and ±3σwith: σi,b=σdata,i,b nmodel,i,b ,(10) where nmodel,i,bis the predicted number of counts in bin b and spectrum iand σdata,i,bis the standard deviation of the data in this bin, √ndata,i,b. To investigate the goodness of the fit we generate pseudoexperiments, 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 L(D| N), i.e. the probability of the experimental data Dgiven the set of parameters  Nof our model. We show in Fig. 10 the result for M1,β/γ ,M2and M1,α. The mean of the distributions of M1,β/γ and M2agree well with the value of the data. For M1,α the result demonstrate a modest incompatibility between the data and the model probably arising from an incomplete modelling of αdetector response or an αmiscalibration. This effect is visible for E >6MeVin Fig. 9bottom 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 M1,α. We detail in Sect. 7.4 the results. The pvalue obtained are p=0.38, p=0.04 and p∼0forM1,β/γ ,M2and M1,α respectively. 7.1 SSD and HSD 2νββ decay mechanisms The transition between the ground states of 100Mo and 100Ru, with spin parity 0+, is realized via virtual βtransitions through 1+states of the intermediate nucleus 100Tc. 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 100Tc (HSD hypothesis) [60]. We have found that the SSD mechanism of 2νββ decay to 100Ru ground state reproduces fairly well the data with a p=0.38, while the HSD model does not, p∼0. Since our data clearly favours SSD over HSD mechanism for 2νββ,we have used the SSD model in our final fit. 7.2 Contaminations derived from the fit 7.2.1 LMO crystal contaminations The M1,α spectrum is populated by α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 α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 αregion is the 210Po peak. This peak is largely described by the Q-value component of the crystal bulk. For the 210Po in order to fit as much as possible the particular shape peak in addition to 10 nm, implantation depths of 1 µm and 1 nm are used in the crystals, and implantation depths of 100 nm and 1 µm are used in the Reflectors. In Fig. 11 the left tail of the 210Po 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 238U and 232Th chains are all below 1µ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. 5the peak at 4.8 MeV contains 234U and 226Ra alpha decays. In this analysis this peak is ascribed to 234U, with a significant uncertainty (reported in Table 4)inthe resulting contamination due to the anticorrelation with the 226Ra contribution. Additionally, in this peak we could have a contribution from the neutron capture in 6Li [31]. Neutrons captured in 6Li produce an alpha particle plus tritium, 6Li(n,α)3H, with a total energy 4.782 MeV. We note also that the level of 228Ra is not constrained by any αpeak. 123 675 Page 16 of 24 Eur. Phys. J. C (2023) 83 :675 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 228Th, 226Ra and 40K 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 3in Sect. 4) Component Bulk Posterior Activity from screening [mBq/kg] [mBq/kg] Reflectorsa 238Uto210Pb 9.2±1.0 Refl. only: 0.17 ±0.05 210Pb <17 232Th to 208Pb <2.3 Refl. only: 0.05 ±0.01 Springs 228Ac <217 228Th to 208Pb 20 ±521±5 226Ra to 210Pb 10 ±311±3 40K 3440+450 −340 3600 ±400 Kapton cables 228Ac <139 228Th to 208Pb <28 15 ±10 226Ra to 210Pb <13 8 ±6 Connectorsb 228Ac <442 228Th to 208Pb <339 82 ±38 226Ra to 210Pb <169 15 ±8 Brass screws 228Ac <24 228Th to 208Pb <18 3.5±0.9 210Bic(3.0±0.3)×104620 ±254 Copper supports 228Ac <0.051 228Th to 208Pb <0.052 0.024 ±0.012 226Ra to 210Pb <0.019 <0.04 60Cod0.47 ±0.02 0.04 57Cod0.029 ±0.005 Cryostat screens 228Ac <0.38 228Th to 208Pb <0.40 0.024 ±0.012 226Ra to 210Pb <0.15 <0.04 PE 1Ke 228Ac <4.40.5±0.2 228Th to 208Pb 2.2+2.1 −1.60.3±0.1 226Ra to 210Pb <2.10.65 ±0.08 Screen 300K 228Ac to 208Pb (203+48 −51)mBq 226Ra to 210Pb (94 ±13)mBq 40K(3200 ±400)mBq Reflectorsa 238Uto234U2.7+1.9 −1.6 234U<9.5 230Th <3.5 226Ra to 210Pb 3.4+1.5 −1.2(1.0±0.4)f/(1.7±0.5)g 210Pb to 206Pb 1034+26 −33 232Th <3.9 228Ra to 228Th <504 228Th to 208Pb 2.6+1.4 −1.5(1.1±0.4)g aReflectors take into account all passive elements directly facing the crystals: reflecting foils, PTFE, bonding wires, heaters bConnectors refer to MillMax connectors plus Kapton connectors cThe 210Bi in the Brass Screws accounts for this contamination in all 10 mK (Cables, Connectors, Springs, Copper supports) and infrastructure sources (cryostat screens and PE 1 K) dCo in Copper supports account for this contamination also in Cryostat screens e1K PE accounts for all sources below the 10 mK stage, e.g., 300 K electronics, dilution unit f 214Bi surface measurement with the BiPo-3 detector gExtrapolation from ICPMS measurement, assuming all contamination on surface 123 Eur. Phys. J. C (2023) 83 :675 Page 17 of 24 675 Fig. 9 Experimental data and background model simultaneous fit reconstruction of the 3 CUPID-Mo data spectra. Two upper panels: M1,β/γ ,β/γ’s spectrum with energy deposits in only one detector. Middle panels: M2, multiplicity 2 events, histogram of the 2 summed energies. Two bottom panels: M1,α , multiplicity 1 events in α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 ±1,±2,and ±3σ 123 675 Page 18 of 24 Eur. Phys. J. C (2023) 83 :675 Fig. 10 Distribution of −ln L(D| N)from the toys for the M1,β/γ (left), M2(middle) and M1,α (right) spectra. The red line shows the −ln L(D| N) of the reference fit for each of the spectra Fig. 11 Experimental M1,α spectrum reconstruction showing the components of the M1,α background model fit. Crystal and Reflector contaminations include bulk and surface. The surface contaminations are modelled with an exponential density profile and λ=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 There is clearly a larger 210Po contribution than the rest of the 238U chain, at the level of 96 µBq/kg, possibly introduced during the purification of the enriched material [61]. There are also traces of 190Pt, caused by the crystal growth in a platinum crucible [62] and we find 40K and 90Sr+90Yatthe level of some hundreds of µBq/kg. We note that 210Pb, 87Rb, 90Sr + 90Y and 40K do not represent a potential background for 0νββ search, as the Qβof these radioisotopes is much lower than the 0νββ 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 µm (i.e., including surface contaminations with λ=10 nm and 10 µ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 αregion. The values are larger than the measured radioactivities of the reflectors themselves, in particular, in 226Ra. 123 Eur. Phys. J. C (2023) 83 :675 Page 19 of 24 675 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 214Bi and 208Tl 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,theConnectors,theBrass 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 208Tl. The Screen 300 K accounts for the residual environmental γ’s and the radon present in the gap between the outermost cryostat screen and the external lead shielding. The 226Ra contamination derived from the fit shown in Table 5can be translated into a radon level concentration resulting in (22 ± 3)mBq/m3, which is in good agreement with measurements of 20 mBq/m3provided by the radon mitigation system in the LSM [34]. Figure 12 shows the breakdown of the components in the fit of M1,β/γ . In the region 0.8–3 MeV the dominant contribution is the 2νββ from 100Mo 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νββ region. 7.3 Background index in the 100Mo 0νββ ROI We use our simultaneous fit to reconstruct the background index in the 0νββ region of interest. We chose to calculate the background index in the region ±15 keV around 3034 keV. This range is much wider than the experimental energy resolution, without including any γlines. We sample directly the full posterior distribution produced by JAGS for each step i in the Markov Chain by computing: bi= 67  j=1 Pois(Nj)wi,j E×Mt.(11) Here biis the background index in the 0νββ region of interest, Njis the integral of the spectrum of MC source jin the ROI, wi,jis the weight for source jin step i,Eis 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 bifor 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: b=2.7+0.7 −0.6×10−3counts/keV/kg/year.(12) or, in terms of the number of moles of isotope, moliso, and energy resolution, EFWHM: B=3.7+0.9 −0.8×10−3counts/EFWHM/moliso/year.(13) This is the lowest background index achieved in a bolometric 0νββ 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νββ 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: 8.1+3.5 −2.5×10−5counts/keV/kg/year.(14) 123 675 Page 20 of 24 Eur. Phys. J. C (2023) 83 :675 Fig. 12 Background sources reconstructing the experimental M1,β/γ spectrum, grouped by source location. In blue, 2νββ 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 As shown in Fig. 14, the posterior probability for pile-up allows us to set an upper limit for its background index, < 1.4×10−3counts/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]). Figure15 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: ln(L(D|( N))) = 3  i=1 Nb(i)  b=1 ln (Poiss(ni,b;fi(Eb; N))) +ln (Poiss(NMC j,i,b;ˆ NMC j,i,b)), (15) where, NMC j,i,bis the number of MC events in bin bof source jin spectra i, and ˆ NMC i,bis 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 Fig. 13 Posterior distribution of background index, showing the mode and the smallest 68.3% c.i., 2.7+0.7 −0.6×10−3counts/keV/kg/year largely unchanged with only a small change in the value of the background index. 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 M1,β/γ and M2. 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 M1,α due to the low statistics in each bin of the data; 123 Eur. Phys. J. C (2023) 83 :675 Page 21 of 24 675 – Fit energy region: our reference fit extends from 100 keV to 4 MeV for M1,β/γ 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 M1,α fit: our reference fit extends from [3000–10000] keV. As described at the beginning of Sect. 7the M1,α 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. 17 is naturally marginalising over the uncertainty on the choice of the background sources. Considering all tests in Table 6as a systematic uncertainty (with 2 keV fixed binning) and adding them in quadraTable 6 Background Index in ROI for different fits. The tests allow to check the stability of the model and assess the systematic uncertainties Fit Background index [10−3cts/keV/kg/year] Reference fit 2.7+0.7 −0.6 M1,β/γ threshold = 200 keV 2.8+0.7 −0.6 1 keV fixed binning for M1,β/γ and M22.5+0.6 −0.5 2 keV fixed binning for M1,β/γ and M22.5+0.7 −0.5 20 keV fixed binning for M1,β/γ and M22.9+0.8 −0.6 No crystal bulk 238Uand232Th chains 2.8+0.7 −0.5 No crystal surface 238Uand232Th chains 2.8+0.7 −0.6 No 10 mK sources 238Uand232Th chains 2.2+0.7 −0.5 MC statistics (nuisance parameter) 2.8+0.7 −0.6 M1,α range = 3000–6360 keV 3.8±0.9 ture, the background index in (3034 ±15) keV results in: b=2.7+0.7 −0.6(stat)+1.1 −0.5(syst)×10−3counts/keV/kg/year, (16) or: B=3.7+0.9 −0.8(stat)+1.5 −0.7 (syst)×10−3counts/EFWHM/moliso/year.(17) 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. Fig. 14 Posterior distributions of background index of the several background sources grouped by source location. Also shown is the full posterior distribution 123 675 Page 22 of 24 Eur. Phys. J. C (2023) 83 :675 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 7.5 Residual alpha background Due to our αparticle rejection, background events in the ROI from 226Ra and 228Th subchains in the bulk and the surface of the crystals generally arise only from energy depositions of βor γparticles. However, 238U, 234U, 230Th, 210Po and 232Th could also produce events in the ROI through energy deposits of αparticles. Even if we apply a light yield cut to remove αbackground, it is still possible that some αevents pass this selection cut. From the background index distribution of the crystals, one can separate the background from β/γ decays from that coming from αdecays, as shown in Fig. 16. One can observe that a non-negligible part of the crystal background index is coming from α’s that passes the light yield cut. This αbackground is coming from surface contamination of the crystals. It corresponds to an α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 M1,β/γ 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. Figure17, right, shows the reconstruction of the crystal component of M1,β/γ 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 furFig. 16 Posterior distribution of background index of the crystal from αand β/γ contamination ther 4 LMOs. We expect that in the case where all the crystals face two LDs, as in CUPID, the αbackground contribution should be negligible. 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 ×year exposure. We have performed a simultaneous fit of three data spectra, M1,β/γ , M2and M1,α, to detailed Monte Carlo simulations. The model is performed in a Bayesian framework with a MCMC 123 Eur. Phys. J. C (2023) 83 :675 Page 23 of 24 675 Fig. 17 Left: Experimental M1,β/γ spectrum (in blue) adding a cut to select crystals that face two LDs (in red). Right: Fit reconstruction of the crystal component from M1,β/γ spectrum (in blue), adding a cut to select crystals that face two LDs (in red) approach. We have shown by a fit to a 56Co 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 Li2100MoO4 crystals is sufficient to reach the goals of the future 0νββ experiment CUPID. The radiopurity levels of 226Ra and 228Th are below 0.5 μ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) ×10−3 counts/EFWHM/moliso/year, the lowest in a bolometric 0νββ decay experiment. The detailing of the background achieved in this work enables promising further studies. We can obtain the 2νββ decay rate of 100Mo with high precision. It also allows for studies on various process which could distort the spectral shape, like Bosonic neutrinos, CP violation or 0νββ with Majoron(s) emission. Acknowledgements This work has been performed in the framework of the CUPID-1 (ANR-21-CE31-0014) and LUMINEU programs, funded by the Agence Nationale de la Recherche (ANR, France). We acknowledge also the support of the P2IO LabEx (ANR-10LABX0038) in the framework “Investissements d’Avenir” (ANR-11IDEX-0003-01 Project “BSM-nu”) managed by ANR, France. The help of the technical staff of the Laboratoire Souterrain de Modane and of the other participant laboratories is gratefully acknowledged. We thank the mechanical workshops of LAL (now IJCLab) for the detector holders fabrication and CEA/SPEC for their valuable contribution in the detector conception. F.A. Danevich, V.V. Kobychev, V.I. Tretyak and M.M. Zarytskyy were supported in part by the National Research Foundation of Ukraine Grant No. 2020.02/0011. A.S. Barabash, S.I. Konovalov, I.M. Makarov, V.N. Shlegel and V.I. Umatov were supported by the Russian Science Foundation under Grant No. 18-12-00003. J. Kotila is supported by Academy of Finland (Grant Nos. 314733, 320062, 345869). Additionally the work is supported by the Istituto Nazionale di Fisica Nucleare (INFN) and by the EU Horizon2020 research and innovation program under the Marie Sklodowska-Curie Grant Agreement No. 754496. This work is also based on support by the US Department of Energy (DOE) Office of Science under Contract Nos. DE-AC02-05CH11231, and by the DOE Office of Science, Office of Nuclear Physics under Contract Nos. DE-FG02-08ER41551, DE-SC0011091. This research used resources of the National Energy Research Scientific Computing Center (NERSC) and the IN2P3 Computing Centre. This work makes use of the Diana data analysis software and the background model based on JAGS, developed by the CUORICINO, CUORE, LUCIFER, and CUPID-0 Collaborations. Russian and Ukrainian scientists have given and give crucial contributions to CUPID-Mo. For this reason, the CUPID-Mo collaboration is particularly sensitive to the current situation in Ukraine. The position of the collaboration leadership on this matter, approved by majority, is expressed at https://cupid-mo.mit.edu/collaboration#statement.Majority of the work described here was completed before February 24, 2022. Data Availability Statement This manuscript has no associated data or the data will not be deposited. [Authors’ comment: This manuscript has associated raw and processed data. The data can be made available by signing an agreement with the CUPID-Mo collaboration.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP3.SCOAP 3supports the goals of the International Year of Basic Sciences for Sustainable Development. 123 675 Page 24 of 24 Eur. Phys. J. C (2023) 83 :675 References 1. M. Agostini, G. Benato, J.A. Detwiler, J. Menéndez, F. Vissani Rev. Mod. Phys. 95(2), 025002 (2023). https://doi.org/10.1103/ RevModPhys.95.025002 2. S. Abe et al., Phys. Rev. Lett. 130(5), 051801 (2023). https://doi. org/10.1103/PhysRevLett.130.051801 3. M. Agostini et al., Phys. Rev. Lett. 125(25), 252502 (2020). https:// doi.org/10.1103/PhysRevLett.125.252502 4. D.Q. Adams et al., Nature 604(7904), 53 (2022). https://doi.org/ 10.1038/s41586-022-04497-4 5. O. Azzolini et al., Phys. Rev. Lett. 129(11), 111801 (2022) 6. C. Augier et al., Eur. Phys. J. C 82(11), 1033 (2022). https://doi. org/10.1140/epjc/s10052-022-10942-5 7. G. Anton et al., Phys. Rev. Lett. 123(11), 161802 (2019) 8. I.J. Arnquist et al., Phys. Rev. Lett. 130(11), 062501 (2023) 9. A. Barabash, Universe 6(10), 159 (2020). https://doi.org/10.3390/ universe6100159 10. E. Armengaud et al., Eur. Phys. J. C 77(11), 785 (2017). https:// doi.org/10.1140/epjc/s10052-017-5343-2 11. E. Armengaud et al., Phys. Rev. Lett. 126(18), 181802 (2021). https://doi.org/10.1103/PhysRevLett.126.181802 12. D.Q. Adams, Prog. Part. Nucl. Phys. 122, 103902 (2021) 13. C. Alduino et al., Eur. Phys. J. C 77, 543 (2017) 14. D.Q. Adams et al., Phys. Rev. Lett. 126(17), 171801 (2021). https:// doi.org/10.1103/PhysRevLett.126.171801 15. C. Alduino et al., Eur. Phys. J. C 77(1), 13 (2017). https://doi.org/ 10.1140/epjc/s10052-016-4498-6 16. D. Poda, Physics 3(3), 473 (2021) 17. CUPID pre-CDR (2019). arXiv:1907.09376 [physics.ins-det] 18. C. Augier et al., Phys. Rev. C 107(2), 025503 (2022). https://doi. org/10.1103/PhysRevC.107.025503 19. F.Simkovic,P.Domin,S.V.Semenov,J.Phys.G27, 2233 (2001). https://doi.org/10.1088/0954-3899/27/11/304 20. P. Domin, S. Kovalenko, F. Simkovic, S.V. Semenov, Nucl. Phys. A 753, 337 (2005). https://doi.org/10.1016/j.nuclphysa.2005.03.003 21. Z.G. Berezhiani, A.Y. Smirnov, J.W.F. Valle, Phys. Lett. B 291,99 (1992). https://doi.org/10.1016/0370-2693(92)90126-O 22. R.N. Mohapatra, E. Takasugi, Phys. Lett. B 211, 192 (1988). https:// doi.org/10.1016/0370-2693(88)90832-5 23. C.P. Burgess, J.M. Cline, Phys. Lett. B 298, 141 (1993). https:// doi.org/10.1016/0370-2693(93)91720-8 24. C.P. Burgess, J.M. Cline, Phys. Rev. D 49, 5925 (1994). https:// doi.org/10.1103/PhysRevD.49.5925 25. P. Bamert, C.P. Burgess, R.N. Mohapatra, Nucl. Phys. B 449,25 (1995). https://doi.org/10.1016/0550-3213(95)00273-U 26. C.D. Carone, Phys. Lett. B 308, 85 (1993). https://doi.org/10.1016/ 0370-2693(93)90605-H 27. R.N. Mohapatra, A. Perez-Lorenzana, C.A. de S Pires, Phys. Lett. B 491, 143 (2000). https://doi.org/10.1016/S0370-2693(00)01031-5 28. A.S. Barabash et al., Nucl. Phys. B 783, 90 (2007) 29. S.A. Ghinescu, O. Nitescu, S. Stoica, Phys. Rev. D 105, 055032 (2022) 30. A. Abada, A. Hernandez-Cabezudo, X. Marcano, J. High Energy Phys. 01, 041 (2019) 31. E. Armengaud et al., J. Instrum. 12(08), P08010 (2017). https:// doi.org/10.1088/1748-0221/12/08/P08010 32. E. Armengaud et al., Astropart. Phys. 47, 1 (2013). https://doi.org/ 10.1016/j.astropartphys.2013.05.004 33. M. L’Hour, Revue archéologique de l’ouest 4(1), 113 (1987) 34. R. Hodák et al., J. Phys. G 46(11), 115105 (2019). https://doi.org/ 10.1088/1361-6471/ab368e 35. B. Schmidt et al., Astropart. Phys. 44, 28 (2013) 36. E. Armengaud et al., Eur. Phys. J. C 80(1), 44 (2020). https://doi. org/10.1140/epjc/s10052-019-7578-6 37. C. Alduino et al., Phys. Rev. C 93(4), 045503 (2016). https://doi. org/10.1103/PhysRevC.93.045503 38. O. Azzolini et al., Eur. Phys. J. C 78(9), 734 (2018). https://doi. org/10.1140/epjc/s10052-018-6202-5 39. K. Alfonso et al., Eur. Phys. J. C 82(9), 810 (2022). https://doi.org/ 10.1140/epjc/s10052-022-10720-3 40. O. Azzolini et al., Eur. Phys. J. C 81(8), 722 (2021). https://doi. org/10.1140/epjc/s10052-021-09476-z 41. A. Barabash et al., J. Instrum. 12(06), P06002 (2017) 42. O. Azzolini et al., Eur. Phys. J. C 79(7), 583 (2019). https://doi. org/10.1140/epjc/s10052-019-7078-8 43. S. Rozov et al. (2010). arXiv preprint. arXiv:1001.4383 44. S. Fiorucci et al., Astropart. Phys. 28, 143 (2007). https://doi.org/ 10.1016/j.astropartphys.2007.05.003 45. R. Lemrani, M. Robinson, V.A. Kudryavtsev, M. De Jesus, G. Gerbier, N.J.C. Spooner, Nucl. Instrum. Methods A 560, 454 (2006). https://doi.org/10.1016/j.nima.2005.12.238 46. R. Gurriaran, Personal communication 47. J. Allison et al., Nucl. Instrum. Methods Phys. Res. A 835, 186 (2016) 48. J. Kotila, F. Iachello, Phys. Rev. C 85(3), 034316 (2012). https:// doi.org/10.1103/PhysRevC.85.034316 49. J. Kotila, Personal communication (2012) 50. O.A. Ponkratenko, V.I. Tretyak, Y.G. Zdesenko, Phys. Atom. Nucl. 63, 1282 (2000). https://doi.org/10.1134/1.855784 51. Geant4 Physics Reference Manual 52. A. Gelman, J.B. Carlin, H.S. Stern, D.B. Dunson, A. Vehtari, D. Rubin, Bayesian Data Analysis, 3rd edn. (Chapman & Hall, London, 2013) 53. M. Plummer, in 3rd International Workshop on Distributed Statistical Computing (DSC 2003) (2003), p. 124 54. M. Plummer, JAGS version 3.3.0 user manual (2012) 55. O.G. Polischuk, MDPI Phys. 3(1), 103 (2021). https://doi.org/10. 3390/physics3010009 56. A. Barabash, Universe 6(10), 159 (2020) 57. A. Armatol et al., Phys. Rev. C 104(1), 015501 (2021) 58. D. Chernyak et al., Eur. Phys. J. C 72, 1989 (2012) 59. D.M. Chernyak, F.A. Danevich, A. Giuliani, M. Mancuso, C. Nones,E.Olivieri,M.Tenconi,V.I.Tretyak,Eur.Phys.J.C74, 2913 (2014). https://doi.org/10.1140/epjc/s10052-014-2913-4 60. R. Arnold et al., Eur. Phys. J. C 79(5), 440 (2019). https://doi.org/ 10.1140/epjc/s10052-019-6948-4 61. E. Armengaud et al., J. Instrum. 10(05), P05007 (2015) 62. V. Grigorieva et al., J. Mater. Sci. Eng. B 7, 63 (2017). https://doi. org/10.17265/2161-6221/2017.3-4.002 63. G.Fantinietal.,J.LowTemp.Phys.209, 1024 (2022). https://doi. org/10.1007/s10909-022-02741-9 123