Background model of the CUPID-0 experiment

Abstract

CUPID-0 is the first large mass array of enriched Zn$^{82}$Se scintillating low temperature calorimeters, operated at LNGS since 2017. During its first scientific runs, CUPID-0 collected an exposure of 9.95 kg year. Thanks to the excellent rejection of $\alpha$ particles, we attained the lowest background ever measured with thermal detectors in the energy region where we search for the signature of $^{82}\hbox {Se}$ neutrinoless double beta decay. In this work we develop a model to reconstruct the CUPID-0 background over the whole energy range of experimental data. We identify the background sources exploiting their distinctive signatures and we assess their extremely low contribution [down to $\sim 10^{-4}$ counts/(keV kg year)] in the region of interest for $^{82}\hbox {Se}$ neutrinoless double beta decay search. This result represents a crucial step towards the comprehension of the background in experiments based on scintillating calorimeters and in next generation projects such as CUPID.

1 Introduction

The postulated neutrinoless double beta decay ($0\nu \beta \beta$) consists of two neutrons of an atomic nucleus simultaneously decaying to two protons and two electrons, without the accompanying emission of electron antineutrinos [1]. If observed, $0\nu \beta \beta$ would provide crucial evidence for lepton number violation and it is one of the most sensitive methods to study neutrino properties such as its nature (Dirac or Majorana) and the absolute value of its mass [2]. The experimental signature of $0\nu \beta \beta$ is a peak at the end of the continuous spectrum produced by the electrons emitted in two-neutrino double beta decay, an allowed, although extremely rare, second order nuclear transition. Detectors with excellent energy resolution, such as low temperature calorimeters (historically also called bolometers), are the best candidates to study this process, being able to disentangle the searched peak from the continuous background. However, the energy resolution is only one of the parameters that concur to determine the sensitivity of an experiment for the search of $0\nu \beta \beta$ decay. Others are the number of isotopes under study, the live time, the detection efficiency and the background rate in the energy region of interest (ROI) [3].

The goal of the background model described in this work is to identify the sources of the CUPID-0 background and evaluate their contribution to the ROI around the $^{82}$Se $0\nu \beta \beta$ Q-value ($2997.9\pm 0.3$ keV [4]). This study is, in particular, fundamental for the design of next generation experiments, because the conventional techniques applied to measure radioactivity in materials are not able to probe levels of contamination as low as those required for future $0\nu \beta \beta$ experiments. Therefore the required information must be extrapolated from current rare event experiments.

In this paper, after introducing the CUPID-0 detector and the data production (Sect. 2), we analyze the experimental spectra in wide energy ranges, from a few hundred keV to $\sim$ 10 MeV, to find signatures of background sources (Sect. 3). The energy spectra produced in the detector by each source are then simulated by means of a Monte Carlo code (Sect. 4). The background model (Sect. 5) is constructed by selecting a representative list of sources whose spectra are combined in a Bayesian fit to the experimental data. Information about contaminant activities available from independent measurements or analyses are included through apposite prior distributions. Finally (Sect. 6), we present the fit results, i.e. the activities obtained for the background sources and their contribution to the $0\nu \beta \beta$ ROI, as well as a discussion of systematic uncertainties.

2 CUPID-0: detector and data production

2.1 The CUPID-0 detector

The detection technique used in CUPID-0 experiment [5] is based on cryogenic scintillating calorimeters. These devices allow a simultaneous detection of energy released as heat and light. We exploit both signals to identify different types of interacting particles. The CUPID-0 detector is a five tower array of 26 ZnSe scintillating crystals, 24 enriched in $^{82}\hbox {Se}$ at 95% level and 2 with $^{82}\hbox {Se}$ natural isotopic abundance. The total detector mass is 10.5 kg of ZnSe, equivalent to 5.17 kg of $^{82}\hbox {Se}$. The crystals are interleaved with germanium light detectors (Ge-LDs), that are used to measure the scintillation signal produced in ZnSe by interacting particles. Both ZnSe crystals and Ge-LDs are held in position by means of PTFE clamps and are thermally coupled to a heat bath at $\sim$ 10 mK by means of a copper structure. In order to increase the light collection, each ZnSe crystal is surrounded by a $\hbox {VIKUITI}^{\mathrm{TM}}$ multi-layers reflecting foil produced by 3M. CUPID-0 is hosted in Hall A of Laboratori Nazionali del Gran Sasso (LNGS), inside the cryostat previously used for the CUORICINO and CUORE-0 experiments [6, 7]. The shielding infrastructure is identical to the CUORE-0 one, with the only difference that the 10 mK thermal shield has not been installed, and that the masses of 50 mK and 600 mK shields have been reduced via electrical discharge machining (EDM). The radio-purity of materials used in CUPID-0 experimental setup has been measured with different techniques [8,9,10,11,12], obtaining the results reported in Table 1.

Table 1 Measurements and limits on contaminations of CUPID-0 detector components. Ge-LD radiopurity is certified by UMICORE company. The contaminations of $\hbox {VIKUITI}^{\mathrm{TM}}$ foils have been measured via Inductively Coupled Plasma Mass Spectrometry (ICPMS) at LNGS [9] and with a BiPo detector [12] (private communication). The limits of the other components are taken from Ref. [10]. Error bars are 1 $\sigma$, limits are 90% CL upper limits

Both ZnSe crystals and Ge-LDs are equipped with a neutron transmutation doped (NTD) Ge thermistor [13], working as temperature–voltage transducer. A P-doped Si Joule heater [14, 15], glued to each device, periodically injects a constant energy reference pulse used to measure gain variations induced by temperature fluctuations. The front-end electronics comprises an amplification stage, a six-pole anti-aliasing active Bessel filter and an 18 bits ADC board [16, 17]. The complete data-stream is digitized with a frequency of 1 kHz (2 kHz) for ZnSe (Ge-LD) and saved on disk in NTuples based on the ROOT software framework [18]. A software derivative trigger with channel dependent threshold is implemented online. When a trigger fires on a crystal, the waveforms of the corresponding Ge-LDs are also flagged as signals. For each event on ZnSe we analyze a window of 5 s (1 s before the trigger and 4 s after it). The analysis window of signals on Ge-LDs is 500 ms long (100 ms before the trigger and 400 ms after it). The samples before the trigger provide the baseline temperature of the detector, while the remaining samples are used to determine the pulse amplitude and shape, for evaluating the deposited energy. More details about the CUPID-0 detector construction and performance can be found in Ref. [9] and references therein.

2.2 Data production

This work is based on data collected with CUPID-0 between June 2017 and December 2018, for a total exposure of 9.95 kg year  (Zn$^{82}$Se). Two of the enriched crystals, not properly working, and the two with natural Se are not included in this analysis.

The collected data are processed offline using a C++ based analysis framework originally developed by the CUORE-0 and CUORE collaborations [19]. The specific analysis tools developed for scintillating calorimeters in the framework of CUPID-0 are presented in Refs. [20, 21].

The aim of the data production sequence is to extract from each triggered waveform the corresponding energy release and interaction time. To improve the signal-to-noise ratio, the data are filtered with a software matched-filter algorithm [22, 23]. The filtered amplitude is then corrected for gain instabilities using the reference pulses periodically injected through Si heaters [15]. The corrected amplitude is converted into energy by fitting a parabolic function with zero intercept to the energy of the most intense peaks produced by a $^{232}\hbox {Th}$ source periodically positioned close to the cryostat external shield [20]. The heat released by $\alpha$ and $\beta /\gamma$ of the same energy is slightly different because of the different energy spent in the light channel. To re-calibrate the $\alpha$ events, we identify the most intense $\alpha$ peaks produced by $^{238}\hbox {U}$ and $^{232}\hbox {Th}$ internal contaminations (see Fig. 2), and convert the amplitude to energy using a parabolic function. Data acquired between two calibrations are grouped in a DataSet and processed together through the analysis chain.

We compute time coincidences between detectors within a ±20 ms window, optimized by studying the time distribution of physical coincident events collected during calibrations. Time coincident events are organized in a multiplet structure, which includes information about the triggered crystals and the total energy released in the detector. Since the total event rate is approximately 50 mHz, the probability of accidental (i.e. causally unrelated) coincidences is relatively small ($\sim 10^{-3}$).

Finally, we implement a series of event selection cuts in order to maximize our sensitivity to physics events [20]. Periods of cryostat instability and malfunction are excluded on a crystal-by-crystal basis. A time veto around each event (4 s before and 4 s after) is applied to remove piled up events. We exploit heater pulses to calculate the trigger efficiency (i.e. the probability that an event is detected and reconstructed at the right energy) and the pile-up cut efficiency [19]. We select particle events by requiring a non-zero light signal simultaneously recorded by Ge-LDs. The efficiency of this cut is evaluated by analyzing time coincident events in two crystals, providing a pure sample of particle events, given the negligible probability of random coincidences [10]. The combined efficiency has a constant value of $\varepsilon _{C}=$(95.7 ± 0.5) % above 150 keV.

2.3 Tagging of $\alpha$ particles

Fig. 1

Shape parameter (SP) of light pulses as a function of particle energy. The red solid line is the mean value of $\alpha$ particles SP, while the red dashed line at $\mu _{\alpha }(E)-3\times \sigma _{\alpha }(E)$ is the boundary used to discriminate $\alpha$ from $\beta /\gamma$ events. The blue dashed line at 2 MeV shows the energy below which the particle identification is not applied

The events generated by $\alpha$ particle interaction are identified relying on the different time development of their light pulses with respect to the ones produced by $\beta$/$\gamma$ interactions [24]. Such different pulse shape is quantified by the shape parameter (SP):

$$\begin{aligned} \text {SP} = \frac{1}{Aw_r}\sqrt{\sum _{i=i_M}^{i_M+\omega _r} (y_i-A \cdot S_i )^2} \end{aligned}$$

(1)

where $y_i$ are the samples of the filtered light pulse, A is its maximum amplitude, and $S_i$ are the samples of the filtered average pulse scaled to unitary amplitude and aligned to $y_i$. The summation starts from the index i$_M$ corresponding to the position of the maximum and runs for w$_r$ points ($\sim$ 50) corresponding to the right width at half maximum of S$_i$. The average pulse is made selecting only $\beta /\gamma$ events in the energy range 1.8−2.64 MeV of $^{232}\hbox {Th}$ calibration measurements, with the method described in Ref. [20]. As a consequence, the SP of $\alpha$ events is much higher than the SP of $\gamma$ events. Figure 1 shows the values of SP as a function of energy in CUPID-0 data. Particle identification is difficult below 2 MeV and, thus, it is exploited only above this energy.

Fig. 2

Top left: CUPID-0 $\mathcal {M}_{1\beta /\gamma }$ spectrum with the following peak labeling: (1) $^{65}$Zn, (2) $^{40}\hbox {K}$, (3) $^{208}$Tl. Top right: $\mathcal {M}_{1\alpha }$ spectrum with the following peak labeling: (1) $^{232}\hbox {Th}$, (2) $^{228}\hbox {Th}$, (3) $^{224}\hbox {Ra}$, (4) $^{212}\hbox {Bi}$, (5) $^{212}\hbox {Bi}$ + $^{212}\hbox {Po}$, (6) $^{238}\hbox {U}$, (7) $^{234}\hbox {U}$ + $^{226}\hbox {Ra}$, (8) $^{230}\hbox {Th}$, (9) $^{222}\hbox {Rn}$, (10) $^{218}\hbox {Po}$, (11) $^{214}$Bi + $^{214}\hbox {Po}$, (12) $^{210}\hbox {Po}$, (13) $^{231}\hbox {Pa}$, (14) $^{211}$Bi, (15) $^{147}$Sm. Bottom left: $\mathcal {M}_{2}$ spectrum. Bottom right: $\varSigma _{2}$ spectrum with the same labels used for $\mathcal {M}_{1\beta /\gamma }$ peaks

To discriminate $\alpha$ from $\beta /\gamma$ events, we calculate the mean and the standard deviation of the $\alpha$ particle SP as a function of energy and we set a boundary at $\text {SP}=\mu _{\alpha }(E)-3\times \sigma _{\alpha }(E)$. The $\mu _{\alpha }(E)$ and $\sigma _{\alpha }(E)$ values of SP are calculated excluding the $\beta /\gamma$ events with SP$<6$ (i.e. the cut used in Ref. [5] to select the $\beta /\gamma$ events). The discrimination boundary adopted in this analysis allows to correctly identify energy depositions by $\alpha$ particles with a probability of 99.9% at all energies greater than 2 MeV. This boundary is higher than the $\mu _{\beta /\gamma }(E)+5\times \sigma _{\beta /\gamma }(E)$ contour of $\beta /\gamma$ SP distribution, thus we select $\beta /\gamma$ events with unitary efficiency. The expected number of wrongly identified $\alpha$ events introduces a negligible contamination in the $\beta /\gamma$ spectrum up to $\sim 5$ MeV.

3 Background analysis

Based on the results from previous experiments [10, 25, 26], the sources of background expected in CUPID-0 are:

• the $2\nu \beta \beta$ decay of $^{82}$Se;

• contaminations of the experimental setup (including the detector itself, the cryostat and the shielding) due to ubiquitous natural radioisotopes of $^{232}$Th, $^{238}$U and $^{235}$U decay chains, and $^{40}$K;

• isotopes produced by cosmogenic activation of detector materials, such as $^{60}$Co and $^{54}\hbox {Mn}$ in copper and $^{65}$Zn in ZnSe;

• cosmic muons, environmental $\gamma$-rays and neutrons.

For a better disentanglement of background sources, we exploit the discrimination of $\alpha$ versus $\beta /\gamma$ events and the detector modular design, that allows to tag events producing simultaneous energy depositions in different ZnSe crystals. We then build the following energy spectra:

• $\mathcal {M}_{1\beta /\gamma }$ is the energy spectrum of $\beta /\gamma$ events that triggered only one bolometer (multiplicity 1 events); this spectrum includes also $\alpha$ events with E<2 MeV that, however, provide a minor contribution;

• $\mathcal {M}_{1\alpha }$ is the energy spectrum of multiplicity 1 events produced by $\alpha$ particle interactions;

• $\mathcal {M}_{2}$ is the energy spectrum of events that simultaneously triggered two bolometers (multiplicity 2 events), built with the energies detected in each crystal;

• $\varSigma _{2}$ is the energy spectrum associated to multiplicity 2 events that contains, for each couple of time coincident events, the total energy released in both crystals.

Events with higher order multiplicity are used to evaluate the contribution of muons, that generate electromagnetic showers triggering several bolometers at the same time.

In Fig. 2 we show the $\mathcal {M}_{1\beta /\gamma }$, $\mathcal {M}_{1\alpha }$, $\mathcal {M}_{2}$ and $\varSigma _{2}$ experimental spectra, with labels on the signatures used to identify the background sources.

The main component of the $\mathcal {M}_{1\beta /\gamma }$ spectrum is the continuum produced by the $2\nu \beta \beta$ decay of $^{82}$Se. The most intense $\gamma$ lines exceeding this continuum are produced by the decays of $^{65}$Zn, $^{40}$K, and $^{208}$Tl and are clearly visible also in the $\varSigma _{2}$ spectrum. In Table 2 we report the counting rates of these $\gamma$ lines and of other smaller peaks attributable to the decays of $^{60}$Co, $^{54}$Mn, $^{228}$Ac, and $^{214}$Bi.

Table 2 $\gamma$-Lines identified in the CUPID-0 data and their counting rate in $\mathcal {M}_{1\beta /\gamma }$ and $\varSigma _{2}$ spectra

The peaks observed in the $\mathcal {M}_{1\alpha }$ spectrum are due to $\alpha$ decays occurring in ZnSe crystals or in the detector components directly facing them. These peaks are produced by isotopes belonging to $^{232}$Th, $^{238}$U and $^{235}$U decay chains, and by $^{147}$Sm, a natural long lived isotope with half-life equal to $1.06 \times 10^{11}$ year [27].

In Table 3, we report the counting rates and the energies of the $\alpha$ peaks. All the main peaks in the $\mathcal {M}_{1\alpha }$ spectrum (except the 5.3 MeV line of $^{210}$Po) are centered at the Q-value of the $\alpha$ decays. This means that the corresponding radioisotopes are located in the bulk or near the surface of ZnSe crystals, because the energy of both $\alpha$ and nuclear recoil is detected. The counts in the peaks have been evaluated by means of Gaussian fits with linear background subtraction. Taking into account that the line shape is not perfectly Gaussian and that some peaks are partially overlapped, the rates in Table 3 are affected by systematic error up to $\sim$ 10%. By analyzing the counting rates of the isotopes in each decay chain, we identify the breaking points of secular equilibrium. In $^{232}$Th chain, the progenitor $^{232}$Th has a lower rate with respect to the daughters isotopes, thus we can infer a breaking point at $^{228}\hbox {Th}$ (that can also be at $^{228}\hbox {Ra}$, since there is no signature for this isotope). In $^{238}$U chain, data interpretation is complicated by the fact that the peak at $\sim$ 4.87 MeV includes counts from both $^{234}\hbox {U}$ and $^{226}\hbox {Ra}$. However, the $^{226}\hbox {Ra}$ activity is constrained to be the same of its daughters $^{222}\hbox {Rn}$ and $^{218}\hbox {Po}$, which have relatively short half-life. Therefore, we can infer that the first part of the chain from $^{238}$U to $^{230}\hbox {Th}$ is at equilibrium and that there are two breaking points: the first at $^{226}\hbox {Ra}$ and the second at $^{210}\hbox {Pb}$.

Table 3 Counting rate of $\alpha$-peaks due to $^{238}$U, $^{235}\hbox {U}$ and $^{232}$Th decay chains. All the isotopes appear only as Q-value peak (Q), except $^{210}\hbox {Po}$ that shows also the $\alpha$-peak. Since the line shape is not perfectly Gaussian and some peaks are partially overlapped, the activity evaluation is affected by a systematic error up to 10%, depending on the method used to evaluate the net number of counts in the $\alpha$ peaks. A few lines are depleted by pile-up effects, or summed to $\beta$ coincident events, therefore their activities cannot be directly evaluated

At energies greater than 7.8 MeV, we observe a continuum spectrum with a double bump shape. The first bump is produced by the $^{214}$Bi–$^{214}\hbox {Po}$ decay sequence in $^{238}$U chain, while the second, starting at 8.9 MeV, is due to the $^{212}\hbox {Bi}$–$^{212}\hbox {Po}$ decay sequence of $^{232}$Th chain. In both Bi–Po sequences, the Po half-life is much shorter than the characteristic rise time of thermal pulses (few ms) produced in bolometers, therefore the energy released by Po $\alpha$ decays sums up with the energy deposited by Bi $\beta$ decays. In data production, these events are tagged and calibrated as $\alpha$, because most of the energy is released by the $\alpha$ particle interaction. As a consequence, the energy of $\beta$ component of these events is underestimated by approximately 20%, due to the different calibration of $\alpha$ and $\beta /\gamma$ energy depositions. The Bi–Po signature is present also in the $\mathcal {M}_{2}$ and $\varSigma _{2}$ spectra, because $\beta$ particles and the associated $\gamma$ can cross the reflecting foils around ZnSe crystals and produce $\mathcal {M}_{2}$ events. On the contrary, since the range of $\alpha$ particles is lower than reflecting foil thickness, the $\mathcal {M}_{2}$ and $\varSigma _{2}$ spectra have a small number of events in the range of $\alpha$-lines from 4 to 7 MeV.

4 Monte Carlo simulations

The background sources identified through data analysis are simulated with a Monte Carlo toolkit, called Arby, based on the Geant4 toolkit [28], version 4.10.02. The radioactive decays from the various background sources can be generated in any volume or surface of the CUPID-0 detector, cryostat and shielding implemented in Arby. The primary and any secondary particles are then propagated through the CUPID-0 geometry using the Livermore physics list. The energy deposited in ZnSe crystals is recorded in the Monte Carlo output together with the time at which the interaction occurred. The fraction of energy released by any particle type is also recorded to allow particle identification. Radioactive decays are implemented using the G4RadioactiveDecay database. The decay chains of $^{232}$Th, $^{238}$U, and $^{235}$U can be simulated completely or in part, to reproduce breaks of secular equilibrium. The $2\nu \beta \beta$ simulation is generated under the single-state dominance hypothesis (SSD) in the framework of the Interacting Boson Model (IBM) [29, 30], while the generation of external muons is described in Ref. [31].

In order to implement the detector response function and data production features in the Monte Carlo data, we reprocess the Arby output with a dedicated code. In particular, to account for detector time resolution, we sum energy depositions that occur in the same crystal within a ±5 ms window. The experimental energy resolution is reproduced by applying a Gaussian smearing function with linearly variable width based on measured FWHM of $\gamma$ and $\alpha$ lines. The energy threshold of each detector is modeled with an error function that interpolates the experimental data of trigger efficiency versus energy. These data are collected in dedicated runs in which the heater is used to generate pulses with variable amplitudes (then converted into particle equivalent energies) and the efficiency is calculated, for each pulse amplitude, as the ratio between triggered and generated pulses. Exactly as done in experimental data production, events depositing energy in different crystals within ± 20 ms window are combined into multiplets and pile-up events (see Sect. 2) in the same crystal are discarded. Finally, exploiting the information about the type of particle depositing energy, we reproduce in the Monte Carlo data the same event selection applied in the experimental data to produce the $\mathcal {M}_{1\alpha }$ and $\mathcal {M}_{1\beta /\gamma }$ spectra. Particularly, we include in the $\mathcal {M}_{1\alpha }$ spectrum (with efficiency $>99.9\%$ at E $>2$ MeV) not only the $\alpha$ events, but also the heterogeneous $\beta /\gamma +\alpha$ events due to Bi–Po decay sequences, that in the experimental data are tagged as $\alpha$ events (see Fig. 1).

Fig. 3

Graphic view of the CUPID-0 experimental setup as modeled in Geant4 with the Arby toolkit. External lead and neutron shield are included in the MC but not represented here

To model the cryostat and its shielding we take as reference the scheme developed for the CUORE-0 background model [10], implementing the geometry changes made in CUPID-0. Concerning particle generation, we group together the components that are made of the same material (and thus share equal contaminant concentration) or that cannot be disentangled as they produce degenerate spectra, given the counting statistics of the experimental data. In Fig. 3 we show the geometry of CUPID-0 cryostat and detector as implemented in Arby. The neutron and modern lead (ExtPb) external shields, even if not represented in the figure, are implemented in MC simulations as well (for a detailed scheme and description of these shields see Ref. [10]). The cryostat components where the background sources are generated are the following:

• the Cryostat External Shields (CryoExt) include the Inner Vacuum Chamber (IVC), the super-insulation layers, the Outer Vacuum Chamber (OVC), and the main bath, whose spectra are degenerate;

• the Cryostat Internal Shields (CryoInt) group the 600 mK and the 50 mK shields, that are made of the same copper;

• the Internal Lead Shield (IntPb) is inserted between the IVC and the 600 mK shield and is made of low background ancient Roman lead.

Fig. 4

Comparison among the spectra of a $^{226}\hbox {Ra}$ source simulated in Crystals or in Reflectors, both in the bulk and on the surface of each element. Distinctive Q-value peaks characterize all crystal spectra and, in case of surface simulations, there are also smaller peaks at $\alpha$ energies or a continuum spectrum of degraded $\alpha$s, depending on contaminant depth. Contaminants uniformly distributed in the 70 $\upmu \hbox {m}$ thickness of Reflectors, or simulated with 10 $\upmu \hbox {m}$ depth parameter, give rise to similar step-like shaped spectra due to energy degraded $\alpha$s. Normalization of MC spectra is chosen to allow easy comparison with the experimental data

The CUPID-0 detector itself, reconstructed with high detail in MC simulations, is made of three main components where the background sources are generated:

• the Holder is the supporting structure for the detectors and is made of a special copper alloy (NOSV copper produced by Aurubis company) suitable for cryogenic use and cleaned according to protocols developed in CUORE [32];

• the Crystals are ZnSe cylinders with heights and positions mirroring the real experimental setup;

• the Reflectors are the foils laterally surrounding the crystals. This component is also used to account for the minor contribution from light detectors (see Table 1), from the amount (< 15 cm$^2$/crystal) of copper surface directly facing the edges of ZnSe crystals, and from the other small parts close to the crystals (PTFE spacers, NTDs, and wires), whose spectra are degenerate with those of reflecting foils.

5 Background model

In the construction of a background model, the crucial step is selecting a representative list of sources for fitting the experimental spectra. The analysis of $\alpha$ and $\gamma$ lines presented in Sect. 3 allows to identify the most relevant sources to be included in the model. Alongside these sources, there are other contaminations not producing prominent signatures and whose location (or even emitting isotope) cannot be determined with certainty. In this case, the use of all possible sources introduces too many degrees of freedom in the fit and produces highly correlated results. To avoid this drawback, that would mask the precision of the results, we identify the so-called reference model, a well balanced set of sources, selected according to the above criteria. We then perform some tests in which the source list is modified to investigate the uncertainties related to the choice of background sources.

In the following subsections, we describe the sources used for background model. We distinguish among internal/near and external sources. In the first category, $\alpha$ radiation is not shielded and we must model bulk and surface contamination separately, since they are characterized by different signatures and can produce different counting rates in the $0\nu \beta \beta$ ROI. Conversely, bulk and surface contaminations of external sources produce degenerate spectra and cannot be disentangled.

5.1 Internal/near sources

The internal/near sources are located in Crystals and in detector components directly facing them, modeled by Reflectors. Natural $\alpha$ radiation cannot cross the thickness of reflecting foils and light detectors surrounding the crystals. As a consequence, the $\mathcal {M}_{1\alpha }$ spectrum is made up of events from internal/near sources only.

The breaks of secular equilibrium identified in Sect. 3 are modeled by producing separate Monte Carlo simulations for the different parts of the decay chains, to leave them free to converge on different normalizations when performing the fit.

In Fig. 4 we provide an insight of some $\mathcal {M}_{1\alpha }$ spectra obtained by simulating the decay sequence from $^{226}\hbox {Ra}$ to $^{210}\hbox {Pb}$ in different positions of Crystals and Reflectors. Bulk simulations are obtained by randomly generating the decays inside a volume. Surface contaminations are simulated with exponential density profiles $e^{-x/\lambda }$ (where $\lambda$ is a changeable depth parameter), used to model not perfectly smooth surfaces and diffusion processes of contaminants.

In the background model, we use three types of simulations for modeling the $\mathcal {M}_{1\alpha }$ events produced by contaminations of Crystals:

• bulk, characterized by sharp peaks at Q-value of $\alpha$ decays;

• very shallow with $\lambda =10$ nm, that in addition to the prominent Q-value peaks, exhibit smaller peaks at $\alpha$ energies due to nuclear recoil escapes;

• deep surface with $\lambda =10$ $\upmu \hbox {m}$ (which is approximately the range of natural $\alpha$ particles in ZnSe), that produce the Q-value peaks over a continuum due to degraded $\alpha$ escapes;

Unlike what was possible to do in CUORE-0 background model [10], it is not straightforward to disentangle surface versus bulk contaminations of Crystals, because $\mathcal {M}_{2}$ events produced by $\alpha$-escapes from crystal surfaces are completely absent in CUPID-0 data due to the interposition of reflecting foils. To compensate for the lack of this signature, we developed a new technique based on the time analysis of consecutive $\alpha$ decays. The ratio between the number of parent and time-correlated daughter events releasing all the decay energy in the same crystal depends on the contaminant location. Indeed, given a parent event at the Q-value, the probability to detect a time-correlated event at the daughter Q-value is nearly one in case of bulk contaminations, whereas it is approximately half in case of surface contaminations, because of $\alpha$-escapes from detector surfaces. In particular, for $^{238}$U chain, we count the number of parent events at 5.59 MeV peak of $^{222}\hbox {Rn}$, followed by 6.12 MeV daughter events produced in the same crystal by $^{218}\hbox {Po}$ decay within 3$\times T_{1/2}$ time window (decay scheme in Eq. 2).

(2)

Similarly, for $^{232}$Th chain, we look for time correlated events generated by $^{224}\hbox {Ra}$ (Q = 5.79 MeV) and $^{220}\hbox {Rn}$ $\alpha$ decays (scheme in Eq. 3). In the latter case, we developed a dedicated tool to recover and count the $^{220}\hbox {Rn}$ events which are rejected by standard analysis cuts due to pile-up with $^{216}\hbox {Po}$ decay ($T_{1/2} = 0.145~\text {s}$).

(3)

By combining the experimental data of these time coincidences with Monte Carlo simulations, we constrain the ratio of surface versus bulk contaminations of the middle ($^{226}\hbox {Ra}$–$^{210}\hbox {Pb}$) and lower ($^{228}\hbox {Ra}$–$^{208}\hbox {Pb}$) parts of $^{238}$U and $^{232}$Th decay chains, respectively. The results of this analysis prove that bulk contaminations have higher activity than surface ones. In particular, we obtain that only $\sim$ 15% ($\sim$ 5%) of the parent events at the Q-value of $^{222}\hbox {Rn}$ ($^{224}\hbox {Ra}$) are produced by surface contaminations. In the background model we exploit this information by setting specific priors to constrain the activity of $^{226}\hbox {Ra}$–$^{210}\hbox {Pb}$ and $^{228}\hbox {Ra}$–$^{208}\hbox {Pb}$ surface contaminations relative to the bulk ones.

The other contamination of Crystals (see Table 4 for the complete list) cannot be constrained with this method and, to avoid getting too much correlated results, are modeled as bulk. This choice does not affect the reconstruction of the background at the $0\nu \beta \beta$ ROI.

The simulations we use for modeling the contaminations of Reflectors are:

• very shallow with $\lambda =10$ nm, producing a sharp peak at $\alpha$ energy;

• bulk or deep surface with $\lambda =10~\upmu \hbox {m}$, both characterized by a continuum spectrum due to degraded $\alpha$ particles.

As shown in Fig. 2 and Table 3, the only $\alpha$ line clearly visible in the $\mathcal {M}_{1\alpha }$ spectrum results from $^{210}\hbox {Po}$ decay. Therefore, the only very shallow contamination of Reflectors included in the reference model is the $^{210}\hbox {Pb}$ one. Given the thickness (70 $\upmu \hbox {m}$) and low density ($0.6~\hbox {g/cm}{^3}$) of reflecting foils, bulk and 10 $\upmu \hbox {m}$ surface contaminations produce nearly degenerate spectra. In the reference model we use the 10 $\upmu \hbox {m}$ surface ones. Based on the measured contaminations of reflector foils reported in Table 1, we do not include the upper part of $^{238}$U decay chain (whose contribution is negligible), and we set a prior on $^{232}$Th chain using the the upper limit on $^{228}\hbox {Ra}$ contamination. Therefore, the only unconstrained deep surface contamination of Reflectors is the lower part of $^{238}$U chain, which is split into $^{226}\hbox {Ra}$–$^{210}\hbox {Pb}$ and $^{210}\hbox {Pb}$–$^{206}\hbox {Pb}$, to allow a break of equilibrium.

Fig. 5

Top: comparison between experimental $\mathcal {M}_{1\beta /\gamma }$ spectrum and fit reconstruction. The lower panel shows the bin by bin ratios between counts in the experimental spectrum over counts in the reconstructed one; the corresponding uncertainties at 1, 2, 3 $\sigma$ are shown as colored bands centered at 1. Bottom: same as top for $\mathcal {M}_{1\alpha }$ spectrum

5.2 External sources

The external sources are contaminations in the holder, in the cryostat and in the shields. These sources produce events in the $\mathcal {M}_{1\beta /\gamma }$, $\mathcal {M}_{2}$, and $\varSigma _{2}$ spectra. The $\gamma$ lines that can be used to identify these sources, besides being few, have limited statistics and, thus, cannot be exploited to directly extract information about the position of contaminations. For this reason, we take as reference the background model of the CUORE-0 experiment, that was operated in the same cryostat as CUPID-0. Given the lower statistics of CUPID-0 data and the higher $2\nu \beta \beta$ rate, we have to apply further approximations with respect to the CUORE-0 model. Particularly, we cannot disentangle $^{232}$Th, $^{238}$U and $^{40}$K contaminations in CryoExt from the ones in ExtPb, because their spectra are degenerate. For the same reason, we merge Holder and CryoInt sources. The only exception is $^{54}\hbox {Mn}$, a cosmogenic-activated isotope with $T_{1/2}=312~\text {days}$, that is mainly located in the most recently produced copper of the Holder structure. $^{60}$Co external sources are simulated only in copper components: CryoInt and CryoExt. We use the result from the CUORE-0 background model to set a prior for $^{60}$Co in CryoExt. Moreover, we do not include $^{40}$K contaminations in IntPb shield, since the CUORE-0 model sets an upper limit for this source. Finally, we simulate a $^{210}\hbox {Pb}$ source in ExtPb, because the bremsstrahlung produced by $^{210}$Bi decay was found to produce a sizable amount of events in CUORE-0 experiment and, thus, is expected to provide a contribution also in CUPID-0.

5.3 Environmental sources

The muon flux, even if strongly suppressed by the Gran Sasso rock overburden, is expected to provide a not negligible contribution to the background in the $0\nu \beta \beta$ ROI. Muons interacting in the detector components can produce several $\gamma$ rays triggering high multiplicity events. We exploit this signature to determine the normalization of the simulated muon spectrum. With this method we determine the contribution from muons within a ± 15% systematic uncertainty, depending on the selection of experimental high multiplicity events used to calculate the normalization factor. The obtained result, which is compatible with measurements performed by other experiments [33], is then used to set a prior for muons in the background model.

The contribution due to environmental neutrons and $\gamma$-rays is negligible, as stated in the CUORE-0 background model [10].

Fig. 6

Top: comparison between experimental $\mathcal {M}_{2}$ spectrum and fit reconstruction. The lower panel shows the bin by bin ratios between counts in the experimental spectrum over counts in the reconstructed one; the corresponding uncertainties at 1, 2, 3 $\sigma$ are shown as colored bands centered at 1. Bottom: same as top for $\varSigma _{2}$ spectrum

6 Results

We perform a simultaneous Bayesian fit of $\mathcal {M}_{1\alpha }$, $\mathcal {M}_{1\beta /\gamma }$, $\mathcal {M}_{2}$ and $\varSigma _{2}$ spectra with a linear combination of 33 background sources, to evaluate their activities. We use the JAGS (Just Another Gibbs Sampler) software [34] to define the Bayesian statistical model and to sample the joint posterior Probability Density Function (PDF) of the fit parameters (i.e. the normalization coefficients of the Monte Carlo spectra). More details about the JAGS-based analysis tool for background model fit can be found in Ref. [10]. We use non-negative uniform priors for all fit parameters, with a few exceptions discussed in Sect. 5.

We choose a variable step size binning to minimize the effect of not ideal detector response which manifests itself in line shapes of complicated modeling, especially in the $\alpha$ region. Practically, we define a binning that does not split the peaks in more than one bin. We set the minimum bin step to 15 keV in $\mathcal {M}_{1\beta /\gamma }$ spectrum and to 25 keV in $\mathcal {M}_{2}$ and $\varSigma _{2}$ spectra. We enlarge the binning in the regions with low density of events, in order to minimize the effects of statistical fluctuations.

The fit range extends from 300 to 5 MeV, and from 2 to 11 MeV for $\mathcal {M}_{1\beta /\gamma }$ and $\mathcal {M}_{1\alpha }$ spectrum, respectively. The multiplicity 2 events used to fill $\mathcal {M}_{2}$ and $\varSigma _{2}$ spectra are selected by requiring that both events in the multiplet have energy above a threshold set to 150 keV. Excluding the low energy events allows to bypass problems related to noise pulses that sometimes can be triggered and that are difficult to be discriminated from low energy physics events.

Fig. 7

Pull distribution obtained from the reference fit residuals of all bins. The Gaussian fit to the pull distribution gives $\mu = 0.00 \pm 0.07$ and $\sigma =1.11 \pm 0.06$, with $\chi ^2$/dof = 18.4/15

Fig. 8

Correlation matrix among the reconstructed activities. The source indexes are from Table 4

We label as reference the fit performed with the binning and energy range described in this section, using the sources of the reference background model. The effect of different choices is investigated in the systematic studies. The results of the reference fit to the experimental data collected with a 9.95 kg year Zn$^{82}$Se exposure are shown in Figs. 5 and 6. In these plots, we show the comparison between the experimental and the fit-reconstructed spectra. The pull distribution, obtained from the fit residuals of all bins, is shown in Fig. 7 and is compatible with a Gaussian with $\mu =0$ and $\sigma =1$.

We analyze the marginal posterior distributions of the fit parameters to evaluate the activities of background sources. Most of the marginal PDFs have a Gaussian shape and we calculate their mean and standard deviation to get the activity and its uncertainty. Differently, when the activity of a source is compatible with zero, we quote a 90% upper limit by integrating the posterior PDF.

From the sampling of the joint posterior PDF, we also extract the correlation matrix between the fit parameters, represented in Fig. 8. As expected, the internal/near sources used to fit the $\mathcal {M}_{1\alpha }$ spectrum are not correlated to the others, while the sources representing the same contaminant in different positions of the cryostat are highly (anti-)correlated.

The activities of the sources used in the reference fit are listed in Table 4. These numbers must be read and interpreted keeping in mind the approximations and the choices made in constructing the background model. Particularly, since the fit is performed on the full statistics collected by all CUPID-0 detectors, the activities evaluated for Crystals and Reflectors are average values of real contaminations, that could be not uniformly distributed. Moreover, the representation of the external sources is extremely simplified and does not aspire to establish with accuracy the activities of sources in cryostat and shields. Despite the above caveats, this method allows to constrain the background sources on their specific signatures in the experimental data, and to extrapolate their contribution to the $0\nu \beta \beta$ ROI on a relative scale, independently of the absolute activity evaluated for each source.

Table 4 List of the sources used to fit the CUPID-0 background data in the reference model. The columns contain (1) the name of the component where the source is located, (2) the corresponding mass (or surface), (3) the contaminant, (4) the source index used to plot the correlation matrix in Fig. 8, and (5) the evaluated activity with the statistical uncertainty (limits are at 90% CL). The source activities listed in this table are the result of a model involving some approximations and the corresponding systematic uncertainties are not included here

Fig. 9

Background sources contributing to the $\mathcal {M}_{1\beta /\gamma }$ reconstruction, grouped by source and component. The shaded area corresponds to the 400 keV energy range from 2.8 to 3.2 MeV ($\hbox {ROI}_{bkg}$) chosen to analyze the background in the region of interest around the $^{82}$Se $0\nu \beta \beta$ Q-value. In this plot, the time veto for the rejection of $^{208}$Tl events is not applied, thus the $\hbox {ROI}_{bkg}$ is dominated by the $\beta /\gamma$-events from $^{232}$Th chain contaminations located in Crystals

6.1 Analysis of $0\nu \beta \beta$ ROI and systematics

To analyze the background in the region of interest around the $^{82}$Se $0\nu \beta \beta$ Q-value, we define a 400 keV interval from 2.8 to 3.2 MeV (hereinafter referred to as $\hbox {ROI}_{bkg}$). The energy range chosen for $\hbox {ROI}_{bkg}$ is much wider with respect to the ($23\pm 0.6$) keV FWHM energy resolution at Q-value [5], in order to include enough experimental counts to be used as a benchmark for background model predictions. In Fig. 9, we show the reconstructed $\mathcal {M}_{1\beta /\gamma }$ spectra of different groups of sources, and their contribution to the $\hbox {ROI}_{bkg}$ counting rate.

The counts predicted by the background model in the $\hbox {ROI}_{bkg}$ are 50.5 ± 1.3. This is perfectly compatible with the 52 counts experimentally observed. Most part of $\hbox {ROI}_{bkg}$ events are due to $^{208}$Tl decays inside the Crystals. This radioisotope belongs to the lower part of $^{232}$Th chain and decays via $\beta /\gamma$ with a Q-value of about 5 MeV. Based on the background model results, the $^{228}\hbox {Ra}$–$^{208}\hbox {Pb}$ source in the bulk of Crystals produces $(8.4\pm 0.3)\times 10^{-3}$ counts/(keV kg year). The same contaminant on the surfaces of Crystals results in a rate of $(7.8\pm 1.2)\times 10^{-4}$ counts/(keV kg year). These rates correspond to a total amount of $\sim$ 37 counts in the $\hbox {ROI}_{bkg}$. Nevertheless, since the $^{208}$Tl half-life is relatively short (3.05 min), the $^{208}$Tl events can be rejected by exploiting the time coincidence with its parent, $^{212}\hbox {Bi}$ [20]. For each event in the $\hbox {ROI}_{bkg}$, we check if it is preceded by a $^{212}\hbox {Bi}$-like $\alpha$ event in the same crystal. By applying a 7$\times T_{1/2}$ window time veto to the simulations of $^{228}\hbox {Ra}$–$^{208}\hbox {Pb}$ sources, the predicted counting rates in the $\hbox {ROI}_{bkg}$ become $(3.4\pm 0.6)\times 10^{-4}$ counts/(keV kg year) and $(3.4\pm 0.5)\times 10^{-4}$ counts/(keV kg year) for bulk and surface contaminations of Crystals, respectively. Hence the model predicts the rejection of 34 ± 1 counts. In the experimental data we observe the rejection of 38 events. We ascribe the higher event rejection in the experimental data to random coincidences in the time veto window.

Nonetheless, after applying the time veto, the background model predicts 16.5 ± 0.8 counts in ROI, which is still well compatible with the 14 measured ones.

In Table 5, we report the counting rates reconstructed in the $\hbox {ROI}_{bkg}$ through the background model for the different sources after applying the time veto to reject $^{208}$Tl events. We sum the contributions from surface contaminations at different depths, and from the different components of cryostat and shields, thus reducing the effect of anti-correlation affecting these sources. As shown in Fig. 9, the spectra of background sources in the $\hbox {ROI}_{bkg}$ exhibit either a flat or a decreasing trend, with no significant peaking structures. Since the $\hbox {ROI}_{bkg}$ is symmetrical around $^{82}$Se $\beta \beta$ Q-value, the counting rate in the $\hbox {ROI}_{bkg}$ is a good approximation of the expected rate in the narrower region where the $0\nu \beta \beta$ signature is searched. This is true for all background components except for the $2\nu \beta \beta$ one, because its spectrum has the endpoint at the Q-value of $^{82}$Se $\beta \beta$ decay. The contribution from $2\nu \beta \beta$ source reported in Table 5 is produced exclusively by events with energy $<2950$ keV, while the expected counting rate from $2\nu \beta \beta$ in a 100 keV range centered at $^{82}$Se Q-value is $<3\times 10^{-6}$ counts/(keV kg year).

In order to study the systematic uncertainties of the background reconstruction in the $\hbox {ROI}_{bkg}$, we perform some fits in which the sources are modeled in a different way with respect to the reference fit. Particularly, we performed the following tests:

• 1. a fit with a reduced list of sources in which we exclude the contaminations evaluated as upper limits in the reference fit;

• 2. a fit with Crystals surface contaminations simulated by setting the depth parameter at 0.1 $\upmu \hbox {m}$ instead of 0.01 $\upmu \hbox {m}$;

• 3. a fit in which the $^{226}\hbox {Ra}$–$^{210}\hbox {Pb}$ contamination in Reflectors is removed from the list of sources;

• 4. a fit in which the Reflectors sources simulated with 10 $\upmu \hbox {m}$ depth parameter are replaced by uniformly distributed contaminations;

• 5. a fit in which we add $^{232}$Th and $^{238}$U contaminations on Holder surfaces ($\lambda =10\mu$m), constrained by priors from CUORE-0 background model [10];

• 6. a fit in which we investigate the effect of $^{232}$Th and $^{238}$U surface contaminations on the 50 mK shield surrounding the CUPID-0 tower;

• 7. three fits in which the source list does not include $^{232}$Th and $^{238}$U contaminations in CryoInt, IntPb, and ExtPb, respectively;

Table 5 Counting rates reconstructed in the $\hbox {ROI}_{bkg}$ (from 2.8 to 3.2 MeV) for the different sources, after applying the time veto for the rejection of $^{208}$Tl events. For each value of counting rate we quote first the statistical uncertainty and then the systematic one. In the left column we report the total counting rates from the different components of the experimental setup, while in the right column we provide their breakdown by source. $^{232}$Th and $^{238}$U refers the chain parts producing background in the $\hbox {ROI}_{bkg}$, i.e $^{228}\hbox {Ra}$–$^{208}\hbox {Pb}$ and $^{226}\hbox {Ra}$–$^{210}\hbox {Pb}$, respectively. The counting rate quoted for Reflectors includes also the contribution from surface contaminations of the Holder. The $^{232}$Th source limit (90% CL) corresponds to the maximum one obtained from the analysis of systematic uncertainties. The contribution from $2\nu \beta \beta$ is produced exclusively by events with energy $<2950$ keV

In all of these tests, we obtain pull distributions compatible with a standard Gaussian. Therefore, we analyze the differences in the $\hbox {ROI}_{bkg}$ counting rates to get an estimate of systematic uncertainties, reported in Table 5. We do not quote a systematic uncertainty for $2\nu \beta \beta$ contribution to the $\hbox {ROI}_{bkg}$, because the results from all tests are within a range much smaller than the statistical uncertainty. Crystals surface contaminations are constrained by the time analysis of consecutive $\alpha$ decays. Their counting rate in the $\hbox {ROI}_{bkg}$ has a maximum variation of $\sim$ 30% when fitting with the reduced list (that does not include 10 $\upmu \hbox {m}$ surface contaminations of Crystals) and when performing the tests number 2 and 3.

The systematic uncertainties affecting the $\hbox {ROI}_{bkg}$ counting rate due to Reflectors and Holder contaminations are investigated through tests number 3, 4, and 5. The bulk/deep surface contaminations in Reflectors produce a continuum of degraded $\alpha$ that allows to obtain a good fit to the $\mathcal {M}_{1\alpha }$ spectrum in the [2–4] MeV range. Since $^{232}$Th in Reflectors is constrained by a prior which makes negligible its contribution, $^{226}\hbox {Ra}$–$^{210}\hbox {Pb}$ and $^{210}\hbox {Pb}$–$^{206}\hbox {Pb}$ are the only reflector sources which are left free to fit this continuum. In fit number 3, we investigate the scenario in which Reflectors are contaminated only by $^{210}\hbox {Pb}$–$^{206}\hbox {Pb}$. The result is that the experimental counts in the [7–7.5] MeV range of $\mathcal {M}_{1\alpha }$ spectrum are not reconstructed by the model and we get a 5$\sigma$ disagreement in that bin. We conclude that a contribution from the $^{226}\hbox {Ra}$–$^{210}\hbox {Pb}$ source in Reflectors is needed to preserve the fit quality and we estimate that its activity (and thus its counting rate in the ROI) must be at least half of that evaluated in the reference fit. On the other hand, the result of fit number 4, in which we choose an equally plausible model for the distribution of contaminants in Reflectors is that the $\hbox {ROI}_{bkg}$ counting rate from this source increases by $\sim$ 50%. The fit number 5 is aimed at investigating the effect of contaminations on Holder surfaces. This background source does not have a specific signature in the experimental data, because most of the $\alpha$ particles generated at Holder surfaces are absorbed by Reflectors. However, some $\beta$ particles from $^{238}$U and $^{232}$Th chains can cross the reflector foils and produce events in the $\hbox {ROI}_{bkg}$. Since the Holder is made of the same NOSV copper used in CUORE-0, we exploit the values of $^{232}$Th and $^{238}$U surface contaminations measured with CUORE-0 detector to constrain these sources. The result is that the reconstructed rate in the $\hbox {ROI}_{bkg}$ increases by $\sim 1.5\times 10^{-4}$  counts/(keV kg year). Particularly, the upper limit on the background due to $^{232}$Th in Reflectors and Holder becomes less stringent: $<3.3\times 10^{-4}$ counts/(keV kg year). Similarly, in test number 6 we evaluate the systematic uncertainty affecting the $\hbox {ROI}_{bkg}$ reconstruction in the case that $^{232}$Th and $^{238}$U surface contaminations on the 50 mK shield surrounding CUPID-0 tower are not negligible with respect to the bulk contaminations of CryoInt. Also in this case, as expected, the fit predicts a higher rate in the $\hbox {ROI}_{bkg}$, with an increase of $\sim 4\times 10^{-4}$  counts/(keV kg year) with respect to the reference one. The fits of test number 7 are used to investigate how the uncertainty on location and description of sources in cryostat and shields is propagated to the estimate of their contribution to the $\hbox {ROI}_{bkg}$.

Finally, we performed some tests in which we varied the minimum bin size (from 5 to 50 keV) and the energy calibration (according to the residuals of a $^{56}\hbox {Co}$ calibration measurement reported in Ref. [20]). The changes in source activities and $\hbox {ROI}_{bkg}$ reconstruction are much smaller than the uncertainties quoted in Tables 4 and 5, thus the corresponding systematic errors are negligible.

7 Conclusion and perspectives

In this paper we fit the CUPID-0 data using 33 radioactive sources, modeled via Monte Carlo simulations. We identify the contribution of the various background sources to the $\hbox {ROI}_{bkg}$ counting rate and we perform an analysis of the corresponding systematic uncertainties. Excluding the $2\nu \beta \beta$ decay contribution (which is negligible at the Q-value of $^{82}$Se $\beta \beta$ decay), we obtain that $\sim$ 44% of background rate in the $\hbox {ROI}_{bkg}$ is produced by cosmic muon showers, while the remaining fraction is due to radioactive decays in Crystals ($\sim$ 33%), in Reflectors and Holder ($\sim$ 6%), and in Cryostat and Shields ($\sim$ 17%).

Based on these results, an upgrade of the CUPID-0 detector has been scheduled in order to reduce the background level in the ROI and to further improve the comprehension of background sources. In CUPID-0 Phase-II we plan to install a muon veto, which will be implemented through a system of plastic scintillators in the external experimental setup. Moreover we will investigate the effect of removing reflecting foils. This will also allow us to get more information about surface contaminations of Crystals, through the analysis of $\mathcal {M}_{2}$ and $\varSigma _{2}$ spectra, as performed in CUORE-0 [11]. Finally, we will add an ultra-pure copper vessel at 10 mK, acting as thermal and radioactive shield, that is expected to further reduce the counting rate due to contaminations of cryostat and shields.

The CUPID-0 results on $\alpha$-background rejection, further strengthened by the analysis presented here, are a founding pillar of the next-generation CUPID experiment [35, 36], based on scintillating calorimeters.