Background identification in cryogenic calorimeters through $\alpha -\alpha$ delayed coincidences

Abstract

Localization and modeling of radioactive contaminations is a challenge that ultra-low background experiments are constantly facing. These are fundamental steps both to extract scientific results and to further reduce the background of the detectors. Here we present an innovative technique based on the analysis of $\alpha -\alpha$ delayed coincidences in ${}^{232}$Th and ${}^{238}$U decay chains, developed to investigate the contaminations of the ZnSe crystals in the CUPID-0 experiment. This method allows to disentangle surface and bulk contaminations of the detectors relying on the different probability to tag delayed coincidences as function of the $\alpha$ decay position.

1 Introduction

Fig. 1

Pictures of the CUPID-0 detector. From left to right: a ZnSe crystal, the same crystal surrounded by the reflecting foil, the Ge light detector mounted on top, the CUPID-0 array of 26 scintillating calorimeters

Experiments searching for rare events, such as neutrinoless double-beta ($0\nu \beta \beta$) decay [1], demand for a detailed background understanding in order to implement possible reduction techniques and to collect crucial information for next-generation experiments [2]. Depending on the detector features, different techniques are adopted to perform this analysis, exploiting all the information contained in the data themselves, such as particle energy, event topology, time correlation, and particle type. Cryogenic calorimeters [3, 4] developed to search for the $0\nu \beta \beta$ decay have already demonstrated how to use these techniques to understand and reduce the background. For example, thanks to the experience gained with the Cuoricino experiment [5], a new conceptual design was introduced for the detector holder, thus mitigating the background due to degraded-energy $\alpha$-particles emitted by copper contaminations in the CUORE-0 experiment [6]. CUORE-0 in turns provided the first background model for cryogenic calorimeters [7], on which the CUORE background budget [8] is based. Over the past decade, we reached a deeper understanding of the background and we further reduced it through scintillating calorimeters [9,10,11,12], which introduced the groundbreaking possibility to identify the interacting particles.

CUPID-0 is the first 10 kg-scale demonstrator of such technique and allowed to reach the lowest background ever measured by cryogenic calorimeters, i.e. $3.5 \times 10^{-3}$ counts/ $(\text {keV\,kg\,yr})$ in the region of interest around the ${}^{82}$Se $\beta \beta$ decay Q value ($Q_{\beta \beta }$= 2997.9 ${\pm }$ 0.3 keV [13]), characterized by an average energy resolution of (20.05 ${\pm }$ 0.34) keV FWHM [14, 15]. Such impressive low background rate was achieved by combining the $\alpha$-particle identification (and rejection) with the analysis of time-delayed coincidences. In particular, we tagged potential ${}^{212}$Bi $\alpha$ decays and we vetoed any event occurring within 7 half-lives of its daughter ${}^{208}$Tl ($T_{1/2}=3.05$ min), that $\beta$ decays with a high Q value (5 MeV). In this way, we reduced the background in the region of interest by a factor $\sim$4, at the cost of only 6% dead time [14].

In general, event-tagging based on time-correlations is a widespread tool to identify, quantify, and reduce the background of rare event experiments [16, 17]. Therefore, we decided to analyze the delayed coincidences due the $\alpha$-decay sequences in ${}^{232}$Th and ${}^{238}$U chains to improve the CUPID-0 background model [18]. Indeed, the $\alpha$-decay features, especially the short range of energy deposition, allow to study the contaminant localization. In this paper, we describe how we analyzed the $\alpha -\alpha$ delayed coincidences in the CUPID-0 data to extract information about the position (bulk vs surface) of crystal contaminations. This is very important to help designing next-generation bolometric experiments searching e.g. for $0\nu \beta \beta$ decay, because the background index induced by such contaminations strongly depends on the their location.

2 Experimental setup

CUPID-0 is the first 10 kg-scale CUPID [2] demonstrator using enriched scintillating calorimeters to search for $0\nu \beta \beta$ decay of ${}^{82}$Se. The CUPID-0 detector is an array of 24 Zn${}^{82}$Se crystals 95% enriched in ${}^{82}$Se and two ZnSe crystals with natural Se, for a total mass of 10.5 kg. When a particle interacts in a ZnSe crystal, it produces a measurable temperature rise proportional to the energy deposit, and a light emission that allows for particle identification. The typical rise and decay times of signal pulses in ZnSe crystals are 14 ms and 36 ms, respectively [19]. The ZnSe crystals are held in a copper frame through small PTFE clamps and laterally surrounded by 70 $\upmu$m thick Vikuiti${}^\mathrm{{TM}}$ reflective foil to enhance light collection. To measure the light signal, 170 $\upmu$m thick germanium wafers operated as calorimetric detectors [20] are faced to the ZnSe crystals. Both light detectors and ZnSe crystals are equipped with a Neutron Transmutation Doped (NTD) Ge thermistor [21], acting as temperature–voltage transducer. The detector is operated at a base temperature of $\sim$10 mK in an Oxford 1000 ${}^3$He/${}^4$He dilution refrigerator located underground in the Hall A of the Laboratori Nazionali del Gran Sasso (Italy). The reader can find some pictures of the detector in Fig. 1 and more details in Ref. [19].

3 Data production

Fig. 2

Search for delayed coincidences in the ${}^{238}$U (left) and ${}^{232}$Th (right) decay chains. The grey spectrum comprises the not piled-up ${\mathcal {M}}_{1}$ $\alpha$ events. We tag as daughter (blue) all the events within a $5 \; T_{1/2}$ time-window after a candidate parent event recorded at the ${}^{222}$Rn (left) or ${}^{224}$Ra (right) Q value peak (red). The spectrum of ${}^{224}$Ra daughters (right, in blue) is miscalibrated due to the ${}^{220}$Rn$-^{216}$Po pile-up. These events are rejected by the pileup rejection cuts, therefore they are not included in the ${\mathcal {M}}_{1\alpha }$ spectrum (grey)

In this work, we analyze the spectrum of $\alpha$-particles detected by CUPID-0 Phase I, which lasted from June 2017 to December 2018 with a live-time of 74% for physics runs. The continuum data stream from ZnSe detectors is amplified and filtered with a 120 dB/decade, six-pole anti-aliasing active Bessel filter and saved on disk with a sampling frequency of 1 kHz by a custom DAQ software package [22]. We run a derivative trigger to identify the heat pulses and save a 5 s window for each detected signal. We save the trigger timestamp of each event with a 1 ms precision, given by the sampling frequency. To get the energy deposited in each event, we extract the pulse amplitude by applying a software matched-filter [23], which improves the signal-to-noise ratio and, thus, the energy resolution. We convert the pulse amplitude into energy by fitting a parabolic function with zero intercept to the energy of the most intense $\alpha$-peaks produced by ${}^{238}$U and ${}^{232}$Th internal contaminations of ZnSe crystals [18]. We select particle events by requiring a non-zero light signal simultaneously recorded by light detectors and we tag the $\alpha$-particles relying on the light pulse shape parameter defined in Ref. [24], which allows to discriminate >99.9% of $\alpha$ events from $\beta /\gamma$ ones at energies >2 MeV. We tag the events that simultaneously trigger more than one crystal within a ${\pm }$20 ms time window, assigning a multiplicity label (${\mathcal {M}}_\#$) equal to the number (#) of crystals hit. Since the total event rate is approximately 50 mHz, the probability of accidental coincidences is almost negligible ($\sim 10^{-3}$). Finally, we analyze the waveform of each triggered event to label piled-up events (1 s before and 4 s after trigger) for which the energy reconstruction is not reliable. The data from two enriched crystals, not properly working [19], and from the two natural crystals are not considered in the current analysis, therefore the total active mass of the detector is 8.74 kg, with a corresponding exposure of 9.95 kg yr.

4 Search for delayed coincidences

In CUPID-0, the $\alpha$-particles are not able to pass through the reflecting foils surrounding the ZnSe detectors, therefore we search for time-correlated events occurring in the same crystal. The factors that mostly affect the capability to correctly identify delayed coincidences are the time resolution of the detector and the event rate (r). The first sets a constraint on the minimum half-life of the daughter nuclide that allows for the two events to be resolved in time. The latter, together with the time window opened to search for delayed coincidences ($\varDelta t_w$), determines the probability of random coincidences:

$$\begin{aligned} P_\text {random} = 1-e^{-r \varDelta t_w} \simeq r \varDelta t_w \quad \quad (\varDelta t_w \ll 1/r) \end{aligned}$$

Since $P_\text {random}$ has to be kept $\ll$1 and $\varDelta t_w$ must be chosen of the order of a few half-lives of the daughter nuclide ($T_{1/2}$), the event rate restricts the possibility of searching for delayed coincidences in isotopes with $T_{1/2} \ll 1/r$ only.

In CUPID-0, the detector time resolution is of the order of few ms and the rate of $\alpha$ events is at maximum $1.7\times 10^{-4}$ Hz/crystal and, on average, $6.3\times 10^{-5}$ Hz/crystal. Therefore, the most suitable $\alpha$-decay sequences in ${}^{238}$U and ${}^{232}$Th chains for this analysis are:

In order to search for delayed coincidences produced by crystal contaminations, we process the data as follows.

• 1. We label as candidate parents ($N_P$) all the single-hit (${\mathcal {M}}_{1}$) not piled-up $\alpha$-events at the Q value peak of the first decay in the sequence, within a ${\pm } 1.5 \, \sigma$ energy resolution range. This is a good compromise to select a large fraction of candidate parents, without including too much background from the continuum underlying the peaks (see red histograms in Fig. 2). This selection focuses the analysis on crystal contaminations, being the only ones that can produce a signal event at the Q value.

• 2. At each candidate parent, we tag as daughters all the events occurring in the same crystal within a time window $\varDelta t_w = 5 \, T_{1/2}$ of the second decay (blue histograms in Fig. 2). The length of the coincidence window is optimized to select a large fraction of signal ($\sim$97%), while keeping random coincidences at a negligible level. We only require that a daughter is an $\alpha$-event, without applying multiplicity and pile-up cuts. In this way, we can identify a delayed coincidence even if an uncorrelated event simultaneously triggers another detector or if pile-up occurs. The latter case is particularly frequent in the ${}^{220}$Rn$-^{216}$Po decay sequence.

• 3. When two candidate parent events occur in the same detector within $\varDelta t_w$, we discard both parents and their daughters from the analysis. In this way, we reduce the contribution from random delayed coincidences and we avoid ambiguity in the assessment of the $\varDelta t$ between couples of parent-daughter events. The expected number of these random coincidences between parent candidates is given by:

  $$\begin{aligned} N_{PP}=\sum _{ch}N_{P}^{ch} P_\text {random} \simeq \sum _{ch} N_{P}^{ch} r_{P}^{ch} \varDelta t_{w} \end{aligned}$$

  (1)

  where $N_{P}^{ch}$ is the number of candidate parent events detected by each channel ch, and $r_{P}^{ch}$ is their rate. In Table 1 (last row), we check that the number of random coincidences between parent events found in the data ($N^{obs}_{PP}$) is compatible with the expected value calculated through Eq. (1), finding a very good agreement in both chains.

The energy spectra of parent and daughter events resulting from this analysis are shown in Fig. 2, together with the spectrum of the ${\mathcal {M}}_{1}$ $\alpha$-events passing the pile-up rejection cut (${\mathcal {M}}_{1\alpha }$).

In the search for ${}^{222}$Rn −${}^{218}$Po delayed coincidences belonging to the ${}^{238}$U chain, the spectrum of daughter events exhibits a clear peak at the ${}^{218}$Po Q value (Fig. 2 (left)), demonstrating the effectiveness of this technique. Since the time window used in this case is relatively long (15.5 min), we also observe random daughter events corresponding to a fraction of $\sim$1% of the ${\mathcal {M}}_{1\alpha }$ spectrum. To determine the number of detected delayed coincidences ($N_C$), i.e. couples of time-correlated parent–daughter events, we exploit the daughter energy signature and we count the number of events falling in a ${\pm } 3 \sigma$ energy resolution range centered at the ${}^{218}$Po Q value.

In the ${}^{232}$Th chain we have a different situation due to the pile-up between ${}^{220}$Rn and ${}^{216}$Po events. Indeed, when searching for daughters of ${}^{224}$Ra decay, most of the selected events are tagged as piled-up and their energies are misreconstructed by the standard data processing (which just discards them). This is why in Fig. 2 (right), where we plot all daughter events including those miscalibrated due to the pile-up, we observe a continuous bump instead of two peaks at the ${}^{220}$Rn and ${}^{216}$Po Q values. Nevertheless, having a good energy reconstruction of these events is not essential, because the information about time correlation is sufficient for the goal of the analysis presented hereafter, which requires $N_C$ to be determined. For this purpose, we simply count the number of ${}^{224}$Ra events followed by a ${}^{220}$Rn −${}^{216}$Po $\alpha$–$\alpha$ piled-up event, that provides an unambiguous signature to identify this sub-chain of 3 consecutive $\alpha$ decays. We conservatively discard the piled-up events spaced less than 80 ms in time because above this threshold we are able to precisely trace back the second pulse amplitude to a full-energy ${}^{216}$Po decay deposition.

Table 1 Summary of the parameters used to identify delayed coincidences and the corresponding numerical results. As discussed in Sect. 5, the ratio between $N_P$ and $N_C$ depends on the contaminant position. Moreover, in the ${}^{232}$Th chain, the daughter selection is further constrained to detect 3 consecutive $\alpha$-decays, the third occurring with a $\varDelta t>80$ ms

In Table 1, we summarize the parameters used to search for delayed coincidences in ${}^{238}$U and ${}^{232}$Th chains and we report the corresponding results obtained for $N_P$ and $N_C$. As observed in previous studies [19, 25], contaminations are not homogeneously distributed over all detectors because of an increasing improvement of their radiopurity in the different production batches. Thus, the number of observed delayed coincidences in the different crystals reflects this inhomogeneity. Nevertheless, we find that the $N_C/N_P$ ratio is nearly constant in almost all crystals.

In order to check our selection of delayed coincidences and quantify the amount of random ones, we analyze the time distribution of the $\varDelta t$ between couples of parent–daughter events. Indeed, the $\varDelta t$ of physical time-correlated events follows an exponential distribution with a characteristic time parameter equal to the mean-life of the daughter, whereas the $\varDelta t$ distribution of random coincidences can be approximated as flat when $\varDelta t_w \ll 1/r$.

As shown in Fig. 3, the measured $\varDelta t$ of the delayed coincidences identified in the ${}^{238}$U (left) and ${}^{232}$Th (right) chains are distributed with an exponential profile compatible with the half-lives of ${}^{218}$Po [26] and ${}^{220}$Rn [27], respectively. The flat background results to be compatible with zero in both cases, pointing us out a negligible number of random coincidences. This is also confirmed by calculating the expected value of random coincidences:

$$\begin{aligned} N_{rnd} \simeq \sum _{ch} (N_P^{ch} - N_{C}^{ch})\, r_{D}^{ch} \varDelta t_{w} \end{aligned}$$

(2)

where $r_{D}^{ch}$ is the $\alpha$-event rate of unpaired daughter-like events (i.e. having the same signature of daughters in term of energy for ${}^{218}$Po or pile-up structure for ${}^{220}$Rn −${}^{216}$Po) not in delayed coincidence with a parent. According to Eq. (2), $N_{rnd} \lesssim 1$ in both ${}^{238}$U and ${}^{232}$Th decay chain analysis.

Fig. 3

Distribution of the $\varDelta t$ between couples of parent–daughter events in the ${}^{238}$U (${}^{222}$Rn–${}^{218}$Po, left) and ${}^{232}$Th (${}^{224}$Ra–${}^{220}$Rn, right) decay chains. The fit function (solid line) is composed by an exponential plus a flat background (dashed line at zero counts) to account for possible random delayed coincidences, whose integral ($\text {N}_{\text {BKG}}$) eventually results to be compatible with zero in both cases. The half-life parameter reconstructed by the fit is compatible in both cases with the values reported in literature, thus confirming the effectiveness and the reliability in the identification of delayed coincidences

5 Localization of crystal contaminations

In the experiments searching for rare events, it is fundamental to localize the radioactive contaminations because the background rejection techniques have different efficiencies depending on their position. In cryogenic calorimeters, surface contaminations are of particular concern because the whole crystal volume is sensitive in detecting particle interaction, without any dead-layer. A very effective way traditionally used to identify surface contaminations of cryogenic calorimeters consists in analyzing the spectrum of ${\mathcal {M}}_{2}$ events comprised of $\alpha$-recoil coincidences in neighbours crystals [7]. This method cannot be applied to CUPID-0 Phase I data analysis, because the reflective foil around the ZnSe crystals absorbs the $\alpha$-particles escaping from their surfaces, preventing the detection of $\alpha$-recoil ${\mathcal {M}}_{2}$ events. To overcome such limitation, we conceived an innovative method based on the analysis of $\alpha$–$\alpha$ delayed coincidences. As shown in the next section, both bulk and surface contaminations produce candidate parent events at the Q value, allowing to search for delayed coincidences with the procedure introduced in Sect. 4. Since the ratio between the number of detected delayed coincidences ($N_C$) to the number of candidate parents ($N_P$) depends on the source location, we can extract information about it. As sketched in Fig. 4, if the contamination is in the crystal bulk, the probability to detect a full-energy daughter event given a candidate parent observed at the Q value, $p \left( D_Q|P_Q \right)$, is almost 1. Conversely, when contaminations are on crystal surfaces, the $\alpha$-particles have a not negligible probability to escape the detector. Thus in this case, the conditional probability $p \left( D_Q|P_Q \right)$ to detect a daughter event at the Q value is significantly $<1$.

Fig. 4

Sketch of $\alpha -\alpha$ delayed coincidences for bulk (left) and surface (right) crystal contaminations. Parent and daughter decays are represented in red and blue, respectively. In the bulk case, there is nearly a 1:1 ratio between detected daughter events and parent candidates. This ratio falls below 1, when contaminations are on the detector surfaces due to the escape of the $\alpha$ emitted in the daughter decay

In order to determine the activity ratio $r = A_s / A_b$ between surface (s) and bulk (b) contaminations of a particular decay sub-chain, we solve the following system of equations:

$$\begin{aligned} {\left\{ \begin{array}{ll} N_P = A_b \, T (\varepsilon _P^b + r \, \varepsilon _P^s)\\ N_C = A_b \, T (\varepsilon _C^b + r \, \varepsilon _C^s) \end{array}\right. } \end{aligned}$$

(3)

in which the number of candidate parents ($N_P$) and the number of delayed coincidences ($N_C$) are expressed as a function of the contaminant activities (A), the measurement livetime (T), and the detection efficiencies ($\varepsilon$). In this system, the different value of $p \left( D_Q|P_Q \right)$ exploited to disentangle bulk and surface contaminations (hereafter labeled as $p_C$) enters in the $\varepsilon _C$ terms, which can be expressed as $\varepsilon _C = \varepsilon _P \, p_C$.

By solving the system in Eq. (3), we finally obtain the formula to calculate r:

$$\begin{aligned} r = \frac{\varepsilon _P^b \left( N_P \, p_C^b - N_C \right) }{\varepsilon _P^s \left( N_C - N_P \, p_C^s \right) } \end{aligned}$$

(4)

The physical constraint to be respected in order to get positive results for r is:

$$\begin{aligned} p_C^s \le N_C / N_P \le p_C^b \end{aligned}$$

(5)

This is consistent with the fact that in an experiment we can expect to observe delayed coincidences from a combination of bulk and surface contaminations. If one of them is dominant, the experimental ratio $N_C / N_P$ will approach the range limits.

6 Evaluation of delayed coincidence probability

Fig. 5

Monte Carlo spectra of ${}^{238}$U decay chain, zoomed on ${}^{222}$Rn and ${}^{218}$Po peaks, obtained by simulating the contaminants in the crystal bulk (left), and on crystal surfaces with depth parameters of 10 $\mu$m (center) and 10 nm (right). The fraction of ${}^{222}$Rn–${}^{218}$Po delayed-coincidences (highlighted in red) over the total events recorded at the Q value peaks decreases as the contaminants are simulated in a thinner surface layer due to $\alpha$-particle escapes. The small peaks appearing in the rightmost plot are due to recoil escapes occurring when contaminants are simulated in a very shallow surface layer ($\lambda =10$ nm)

In the previous section, we showed that the ratio between the activity of surface and bulk crystal contaminations can be determined from the experimental data once the efficiencies ($\varepsilon _P^b$ and $\varepsilon _P^s$) and the probabilities to detect delayed coincidences ($p_C^b$ and $p_C^s$) are known. We evaluate these parameters through Monte Carlo simulations.

We simulate the background sources identified in Sect. 4 with a Monte Carlo tool, named Arby, based on the Geant4 toolkit [28], version 10.02. The particles generated by the radioactive decays of interest are propagated using the G4EmLivermore physics list. The decay chains of ${}^{232}$Th and ${}^{238}$U can be simulated completely or in part, to reproduce secular equilibrium breaks. For each energy deposit in the detector, we record the information about the crystal where the interaction took place, the amount of deposited energy, and the time elapsed since the previous event in the decay chain. In order to make the simulation output as similar as possible to the experimental data, we implement the detector response and the data production features in a second step. We associate to each decay an absolute time randomly sampled on the time scale of the experimental data taking, with the exception of the decays occurring within 1 hour from their predecessors in order to preserve the time correlations of interest for this analysis. We account for the detector time resolution by summing up the energy depositions that occur in the same crystal within a ${\pm }5$ ms window, and we group into multiplets the events involving different crystals within ${\pm }20$ ms.

We simulate bulk contaminations by generating the decays in random positions uniformly distributed within the whole crystal volume. The usual approach for simulating surface contaminations in cryogenic calorimeters is to sample the decay positions from an exponential distribution $e^{-x/\lambda }$, where $\lambda$ is the so-called depth parameter. This parameter is usually chosen in the range from few nm up to tens of $\mu$m, in order to reproduce the signatures from shallower to deeper contaminations observed in the experimental data [7, 18]. Different values of the depth parameter can be traced back to different contamination mechanisms. For example, the exposure of a material to the air is expected to produce a very shallow contamination, whereas the treatment of crystal surfaces [29] can originate deeper contaminations.

In Fig. 5 we show the result obtained for the ${}^{238}$U chain simulated in the bulk of crystals and on their surfaces with two different depth parameters. We set the two $\lambda$ values equal to 10 nm and 10 $\mu$m, being much lower and on the same scale of the $\alpha$-particle range, respectively. Even for very shallow surface contaminations, we get a significant fraction of events reconstructed at the $\alpha$-decay Q value. We process the Monte Carlo outputs with the same procedure used to tag the delayed coincidences in the experimental data and we highlight the parent-daughter coincident events in the plot. As expected, the ratio of delayed coincidences over the number of parent candidates decreases as the contamination is simulated in a shallower layer near the crystal surfaces.

The analysis method presented in Sect. 5, provides a single parameter to quantify the activity of surface contaminations, thus a unique model must be chosen to describe them. According to Fig. 5, the deeper surface contaminations (10 $\mu$m) produce a delayed coincidence signal which can be viewed as a combination of a bulk contamination and a shallower surface one. Therefore, in our analysis, we choose the simulations with $\lambda =10$ nm to model surface contaminations and to study the ratio between bulk and surface activities.

Table 2 Detection efficiencies of candidate parents ($\varepsilon _P$) and probabilities of delayed coincidences ($p_C$) evaluated from Monte Carlo simulations of bulk (b) and surface (s) crystal contaminations. For the surface contaminations, we sample the decay positions from an exponential distribution with $\lambda = 10$ nm. Uncertainties are negligible due to the high Monte Carlo statistics

In Table 2, we report the values of $\varepsilon _P^b$, $\varepsilon _P^s$, $p_C^b$, and $p_C^s$ obtained from the MC simulations of ${}^{238}$U and ${}^{232}$Th decay chains. The $\varepsilon _P$ efficiencies are computed by taking into account that a ${\pm } 1.5 \, \sigma$ range was used to select the candidate parents in the experimental data. As expected from a simple geometric reasoning about $\alpha$ escape process, the efficiencies and the probabilities related to surface contaminations are about half with respect to the bulk ones. The only exception is the probability to detect a delayed ${}^{220}$Rn–${}^{216}$Po piled-up event. This is because we are searching for a triple $\alpha$-decay sequence and we have to discard a fraction of piled-up events with $\varDelta t<80$ ms to be consistent with the experimental data processing.

7 Results and discussion

Table 3 Experimental input and final results of the delayed coincidence analysis. Both for ${}^{238}$U and ${}^{232}$Th decay chain we quote: the number of parent candidates ($N_P$), net of background subtraction; the number of detected coincidences ($N_C$) with their binomial uncertainties; the ratio between surface and bulk contamination activities (r) and their absolute values ($A_b$, $A_s$). We quote the r result for ${}^{232}$Th chain with an asymmetric uncertainty range to exclude negative non-physical values

In this section we report the results of the delayed coincidence analysis based on the CUPID-0 data, obtained by combining the experimental data ($N_P$ and $N_C$) with the Monte Carlo evaluations summarized in Table 2. The $N_P$ values reported in Table 1 can not be directly used to calculate the activity ratio r, because they include a fraction of background events. Indeed, other radioactive sources can produce some events falling in the energy range of candidate parents. We exploit the CUPID-0 background model [18] to assess such contribution, which results on the percent scale for both parent peaks in the ${}^{232}$Th and ${}^{238}$U chains. The $N_P$ values obtained after subtracting the expected background counts are reported in Table 3, with an uncertainty that takes into account the Poisson fluctuations. The uncertainty associated to $N_C$ is the Binomial one with $N_C$ successes given $N_P$ trials. After calculating r with Eq. (4), we solve the system in Eq. (3) to get $A_s$ and $A_b$.

The results of this analysis prove that most of CUPID-0 crystal contaminants are located in their bulk. For the ${}^{238}$U sub-chain we get that $\sim$20% of decays occur near crystal surfaces, whereas for the ${}^{232}$Th sub-chain this fraction is constrained in a range between zero and a 13% 1 $\sigma$ upper limit. This information was used to set prior constraints in the CUPID-0 background model [18], allowing for the disentanglement of surface vs bulk crystal contaminations.

Given the total mass ($m=8.74$ kg) and surface ($S=2149$ cm${}^{2}$) of ZnSe crystals used for this analysis, we calculate the specific activities of the $\alpha$-decay sequences for the ${}^{238}$U sub-chain:

$$\begin{aligned} A_b / m&= (16.7 {\pm } 0.5) \; \mu \text {Bq/kg} \\ A_s / S&= (16 {\pm } 4) \; \text {nBq/cm}^2 \end{aligned}$$

and for the ${}^{232}$Th one:

$$\begin{aligned} A_b / m&= (12.4 {\pm } 0.6) \; \mu \text {Bq/kg} \\ A_s / S&= (1.4 ^{+ 5}_{-1.4}) \; \text {nBq/cm}^2 \end{aligned}$$

It is worth noting that, because of secular equilibrium break, these results refer only to the second parts of ${}^{238}$U and ${}^{232}$Th decay chains, which are characterized by higher activities with respect to the first parts (see [18] for more details).

7.1 Systematics discussion

The results of this analysis depend on the efficiencies and probabilities related to surface contaminations reported in Table 2. These parameters are affected by the escape probabilities of $\alpha$s and nuclear recoils. The $\alpha$ escape probability is significantly affected when the contamination depth is changed from $\lambda =10$ nm to a value on the same scale of the $\alpha$ range. For example, by analyzing our data with $\lambda =10$ $\mu$m, the activity ratios r would scale up by a factor $\sim 2$. This confirms that, for our analysis, a deep surface contamination is equivalent to a combination of a bulk and a shallow surface contamination. On the other hand, according to our MC simulations, we can consider the nuclear recoil escape as a second order effect for $\lambda > rsim 10$ nm. Since in the experimental spectrum there are no emerging peaks at the $\alpha$ energies of the isotopes analyzed in this work, we can exclude that shallower contaminations ($\lambda \ll$10 nm) can significantly affect our results.

8 Conclusions

In this work, we presented an innovative analysis technique to study the background sources in cryogenic calorimeters relying on the time-correlation of $\alpha$-decay sequences in ${}^{238}$U and ${}^{232}$Th chains. This method allowed us to disentangle surface and bulk contaminations of ZnSe crystals exploiting the different probability to detect delayed coincidences depending on the contamination depth (see Fig. 5). In particular, we demonstrated that the ${}^{238}$U and ${}^{232}$Th contaminants of CUPID-0 detectors are mainly located in the bulk of crystals. This technique, that was applied for the first time to set prior constraints in the CUPID-0 background model [18], can be adopted also in other experiments for broader purposes. For example, in the analysis of CUORE data [30], delayed coincidences could help to better constrain the background sources [31] and to reject the time-correlated events falling in the region of interest. Moreover, the R&D activities for CUPID [32, 33], CUPID-Mo [34,35,36], and in general the experiments searching for rare events can profit from this technique to study the radioactive contaminations of detector components and to select ultra-pure materials, with also the possibility to analyze other decay sequences in ${}^{238}$U, ${}^{232}$Th and ${}^{235}$U chains [37]. Finally, the analysis of delayed coincidences in CUPID-0 allowed to measure the half-life of ${}^{216}$Po [38].