Background Model of the CUPID-0 Experiment

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 yr. 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}$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 $\sim10^{-4}$ counts/(keV kg yr)) in the region of interest for $^{82}$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.


Introduction
The postulated neutrinoless double beta decay (0νββ) 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νββ 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νββ 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νββ decay. Others are the number of 0123456789().: V,-vol 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νββ Q-value (2997.9 ± 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νββ 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 ∼ 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νββ ROI, as well as a discussion of systematic uncertainties.

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 Se at 95% level and 2 with 82 Se natural isotopic abundance. The total detector mass is 10.5 kg of ZnSe, equivalent to 5.17 kg of 82 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 ∼ 10 mK by means of a copper structure. In order to increase the light collection, each ZnSe crystal is surrounded by a VIKUITI 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 exper- Table 1 Measurements and limits on contaminations of CUPID-0 detector components. Ge-LD radiopurity is certified by UMICORE company. The contaminations of VIKUITI 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 σ , limits are 90% CL upper limits Component 232 Th (Bq/kg) 238 U (Bq/kg) Ge-LD < 6 × 10 −6 < 1.9 × 10 −5 a Limit on 212 Bi-212 Po contamination iments [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.
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 antialiasing 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.

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 Th source periodically positioned close to the cryostat external shield [20]. The heat released by α and β/γ of the same energy is slightly different because of the different energy spent in the light channel. To re-calibrate the α events, we identify the most intense α peaks produced by 238 U and 232 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 (∼ 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 is the boundary used to discriminate α from β/γ events. The blue dashed line at 2 MeV shows the energy below which the particle identification is not applied combined efficiency has a constant value of ε C =(95.7 ± 0.5) % above 150 keV.

Tagging of α particles
The events generated by α particle interaction are identified relying on the different time development of their light pulses with respect to the ones produced by β/γ interactions [24]. Such different pulse shape is quantified by the shape parameter (SP): 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 (∼ 50) corresponding to the right width at half maximum of S i . The average pulse is made selecting only β/γ events in the energy range 1.8−2.64 MeV of 232 Th calibration measurements, with the method described in Ref. [20]. As a consequence, the SP of α events is much higher than the SP of γ 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.   To discriminate α from β/γ events, we calculate the mean and the standard deviation of the α particle SP as a function of energy and we set a boundary at SP = μ α (E)−3×σ α (E). The μ α (E) and σ α (E) values of SP are calculated excluding the β/γ events with SP< 6 (i.e. the cut used in Ref. [5] to select the β/γ events). The discrimination boundary adopted in this analysis allows to correctly identify energy depositions by α particles with a probability of 99.9% at all energies greater than 2 MeV. This boundary is higher than the μ β/γ (E) + 5 × σ β/γ (E) contour of β/γ SP distribution, thus we select β/γ events with unitary efficiency. The expected number of wrongly identified α events introduces a negligible contamination in the β/γ spectrum up to ∼ 5 MeV.

Background analysis
Based on the results from previous experiments [10,25,26], the sources of background expected in CUPID-0 are: -the 2νββ 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 Mn in copper and 65 Zn in ZnSe; -cosmic muons, environmental γ -rays and neutrons.
For a better disentanglement of background sources, we exploit the discrimination of α versus β/γ 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: -M 1β/γ is the energy spectrum of β/γ events that triggered only one bolometer (multiplicity 1 events); this spectrum includes also α events with E<2 MeV that, however, provide a minor contribution; -M 1α is the energy spectrum of multiplicity 1 events produced by α particle interactions; -M 2 is the energy spectrum of events that simultaneously triggered two bolometers (multiplicity 2 events), built with the energies detected in each crystal; -Σ 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 M 1β/γ , M 1α , M 2 and Σ 2 experimental spectra, with labels on the signatures used to identify the background sources.
The main component of the M 1β/γ spectrum is the continuum produced by the 2νββ decay of 82 Se. The most intense γ lines exceeding this continuum are produced by the decays of 65 Zn, 40 K, and 208 Tl and are clearly visible also in the Σ 2 spectrum. In Table 2 we report the counting rates of these γ lines and of other smaller peaks attributable to the decays of 60 Co, 54 Mn, 228 Ac, and 214 Bi.
The peaks observed in the M 1α spectrum are due to α 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 × 10 11 year [27].
In Table 3, we report the counting rates and the energies of the α peaks. All the main peaks in the M 1α spectrum (except the 5.3 MeV line of 210 Po) are centered at the Q-value of the α decays. This means that the corresponding radioisotopes are located in the bulk or near the surface of ZnSe crystals, because the energy of both α 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 ∼ 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 Th (that can also be at 228 Ra, since there is no signature for this isotope). In 238 U chain, data interpretation is complicated by the fact that Table 3 Counting rate of α-peaks due to 238 U, 235 U and 232 Th decay chains. All the isotopes appear only as Q-value peak (Q), except 210 Po that shows also the α-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 α peaks. A few lines are depleted by pile-up effects, or summed to β 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 Po decay sequence in 238 U chain, while the second, starting at 8.9 MeV, is due to the 212 Bi-212 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 α decays sums up with the energy deposited by Bi β decays. In data production, these events are tagged and calibrated as α, because most of the energy is released by the α particle interaction. As a consequence, the energy of β component of these events is underestimated by approximately 20%, due to the different calibration of α and β/γ energy depositions. The Bi-Po signature is present also in the M 2 and Σ 2 spectra, because β particles and the associated γ can cross the reflecting foils around ZnSe crystals and produce M 2 events. On the con-trary, since the range of α particles is lower than reflecting foil thickness, the M 2 and Σ 2 spectra have a small number of events in the range of α-lines from 4 to 7 MeV.

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νββ 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 γ and α 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 M 1α and M 1β/γ spectra. Particularly, we include in the M 1α spectrum (with efficiency > 99.9% at E > 2 MeV) not only the α events, but also the heterogeneous β/γ + α events due to Bi-Po decay sequences, that in the experimental data are tagged as α events (see Fig. 1). 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.
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]; -theCrystals are ZnSe cylinders with heights and positions mirroring the real experimental setup; 2000 3000 4000 5000 6000 7000 8000 9000 10000 11000 Energy (keV)  Fig. 4 Comparison among the spectra of a 226 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 α energies or a continuum spectrum of degraded αs, depending on contaminant depth. Contaminants uniformly distributed in the 70 µm thickness of Reflectors, or simulated with 10 µm depth parameter, give rise to similar step-like shaped spectra due to energy degraded αs. Normalization of MC spectra is chosen to allow easy comparison with the experimental data -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.

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 α and γ 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, α 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νββ ROI. Conversely, bulk and surface contaminations of external sources produce degenerate spectra and cannot be disentangled.

Internal/near sources
The internal/near sources are located in Crystals and in detector components directly facing them, modeled by Reflectors. Natural α radiation cannot cross the thickness of reflecting foils and light detectors surrounding the crystals. As a consequence, the M 1α 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 M 1α spectra obtained by simulating the decay sequence from 226 Ra to 210 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/λ (where λ 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 M 1α events produced by contaminations of Crystals: Eur. Phys. J. C (2019) 79:583 bulk, characterized by sharp peaks at Q-value of α decays; very shallow with λ = 10 nm, that in addition to the prominent Q-value peaks, exhibit smaller peaks at α energies due to nuclear recoil escapes; deep surface with λ = 10 µm (which is approximately the range of natural α particles in ZnSe), that produce the Q-value peaks over a continuum due to degraded α 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 M 2 events produced by α-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 α 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 Qvalue is nearly one in case of bulk contaminations, whereas it is approximately half in case of surface contaminations, because of α-escapes from detector surfaces. In particular, for 238 U chain, we count the number of parent events at 5.59 MeV peak of 222 Rn, followed by 6.12 MeV daughter events produced in the same crystal by 218 Po decay within 3×T 1/2 time window (decay scheme in Eq. 2).
Similarly, for 232 Th chain, we look for time correlated events generated by 224 Ra (Q = 5.79 MeV) and 220 Rn α decays (scheme in Eq. 3). In the latter case, we developed a dedicated tool to recover and count the 220 Rn events which are rejected by standard analysis cuts due to pile-up with 216 Po decay (T 1/2 = 0.145 s). 224 Ra 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 Ra-210 Pb) and lower ( 228 Ra-208 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 ∼ 15% (∼ 5%) of the parent events at the Q-value of 222 Rn ( 224 Ra) are produced by surface contaminations. In the background model we exploit this information by setting specific priors to con-strain the activity of 226 Ra-210 Pb and 228 Ra-208 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νββ ROI.
The simulations we use for modeling the contaminations of Reflectors are: very shallow with λ = 10 nm, producing a sharp peak at α energy; bulk or deep surface with λ = 10 µm, both characterized by a continuum spectrum due to degraded α particles.
As shown in Fig. 2 and Table 3, the only α line clearly visible in the M 1α spectrum results from 210 Po decay. Therefore, the only very shallow contamination of Reflectors included in the reference model is the 210 Pb one. Given the thickness (70 µm) and low density (0.6 g/cm 3 ) of reflecting foils, bulk and 10 µm surface contaminations produce nearly degenerate spectra. In the reference model we use the 10 µ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 Ra contamination. Therefore, the only unconstrained deep surface contamination of Reflectors is the lower part of 238 U chain, which is split into 226 Ra-210 Pb and 210 Pb-206 Pb, to allow a break of equilibrium.

External sources
The external sources are contaminations in the holder, in the cryostat and in the shields. These sources produce events in the M 1β/γ , M 2 , and Σ 2 spectra. The γ 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νββ 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 Mn, a cosmogenicactivated isotope with T 1/2 = 312 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 Cry- oExt. 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 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.

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νββ ROI. Muons interacting in the detector components can produce several γ rays triggering high multiplicity events. We exploit this signature to determine the normalization of the simulated muon spec-trum. 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 γrays is negligible, as stated in the CUORE-0 background model [10].

Results
We perform a simultaneous Bayesian fit of M 1α , M 1β/γ , M 2 and Σ 2 spectra with a linear combination of 33 back-Eur. Phys. J. C (2019) 79:583 ground 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 nonnegative 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 α 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 M 1β/γ spectrum and to 25 keV in M 2 and Σ 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 M 1β/γ and M 1α spectrum, respectively. The multiplicity 2 events used to fill M 2 and Σ 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.
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   Table 4 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 μ = 0 and σ = 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 M 1α spectrum are not correlated to the oth-ers, 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νββ ROI on a relative scale, independently of the absolute activity evaluated for each source.

Analysis of 0νββ ROI and systematics
To analyze the background in the region of interest around the 82 Se 0νββ Q-value, we define a 400 keV interval from 2.8 to 3.2 MeV (hereinafter referred to as ROI bkg ). The energy range chosen for ROI bkg is much wider with respect to the (23 ± 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 M 1β/γ spectra of different groups of sources, and their contribution to the ROI bkg counting rate.
The counts predicted by the background model in the ROI bkg are 50.5 ± 1.3. This is perfectly compatible with the 52 counts experimentally observed. Most part of 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 β/γ with a Q-value of about 5 MeV. Based on the background model results, the 228 Ra-208 Pb source in the bulk of Crystals produces (8.4±0.3)×10 −3 counts/(keV kg year). The same contaminant on the surfaces of Crystals results in a rate of (7.8 ± 1.2) × 10 −4 counts/(keV kg year). These rates correspond to a total amount of ∼ 37 counts in the 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 Bi [20]. For each event in the ROI bkg , we check if it is preceded by a 212 Bi-like α event in the same crystal. By applying a 7×T 1/2 window time veto to the simulations of 228 Ra-208 Pb sources, the predicted counting rates in the ROI bkg become (3.4±0.6)×10 −4 counts/(keV kg year) and (3.4±0.5)×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 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  In this plot, the time veto for the rejection of 208 Tl events is not applied, thus the ROI bkg is dominated by the β/γ -events from 232 Th chain contaminations located in Crystals the expected rate in the narrower region where the 0νββ signature is searched. This is true for all background components except for the 2νββ one, because its spectrum has the endpoint at the Q-value of 82 Se ββ decay. The contribution from 2νββ source reported in Table 5 is produced exclusively by events with energy < 2950 keV, while the expected counting rate from 2νββ in a 100 keV range centered at 82 Se Q-value is < 3 × 10 −6 counts/(keV kg year). In order to study the systematic uncertainties of the background reconstruction in the 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: In all of these tests, we obtain pull distributions compatible with a standard Gaussian. Therefore, we analyze the differences in the ROI bkg counting rates to get an estimate of systematic uncertainties, reported in Table 5. We do not quote a systematic uncertainty for 2νββ contribution to the 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 α decays. Their counting rate in the ROI bkg has a maximum variation of ∼ 30% when fitting with the reduced list (that does not include 10 µm surface contaminations of Crystals) and when performing the tests number 2 and 3.
The systematic uncertainties affecting the 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 α that allows to obtain a good fit to the M 1α spectrum in the [2][3][4] MeV range. Since 232 Th in Reflectors is constrained by a prior which makes negligible its contribution, 226 Ra-210 Pb and 210 Pb-206 Pb are the only reflector sources which are left free to fit this continuum. In fit number Table 5 Counting rates reconstructed in the 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 pro-ducing background in the ROI bkg , i.e 228 Ra-208 Pb and 226 Ra-210 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νββ is produced exclusively by events with energy < 2950 keV 3, we investigate the scenario in which Reflectors are contaminated only by 210 Pb-206 Pb. The result is that the experimental counts in the [7-7.5] MeV range of M 1α spectrum are not reconstructed by the model and we get a 5σ disagreement in that bin. We conclude that a contribution from the 226 Ra-210 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 ROI bkg counting rate from this source increases by ∼ 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 α particles generated at Holder surfaces are absorbed by Reflectors. However, some β particles from 238 U and 232 Th chains can cross the reflector foils and produce events in the 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 ROI bkg increases by ∼ 1.5 × 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 × 10 −4 counts/(keV kg year). Similarly, in test number 6 we evaluate the systematic uncertainty affecting the 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 negli-gible with respect to the bulk contaminations of CryoInt. Also in this case, as expected, the fit predicts a higher rate in the ROI bkg , with an increase of ∼ 4×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 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 Co calibration measurement reported in Ref. [20]). The changes in source activities and ROI bkg reconstruction are much smaller than the uncertainties quoted in Tables 4 and 5, thus the corresponding systematic errors are negligible.

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 ROI bkg counting rate and we perform an analysis of the corresponding systematic uncertainties. Excluding the 2νββ decay contribution (which is negligible at the Q-value of 82 Se ββ decay), we obtain that ∼ 44% of background rate in the ROI bkg is produced by cosmic muon showers, while the remaining fraction is due to radioactive decays in Crystals (∼ 33%), in Reflectors and Holder (∼ 6%), and in Cryostat and Shields (∼ 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 M 2 and Σ 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 α-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.