Background identification in cryogenic calorimeters through \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha -\alpha $$\end{document}α-α delayed coincidences

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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha -\alpha $$\end{document}α-α delayed coincidences in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${}^{232}$$\end{document}232Th and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${}^{238}$$\end{document}238U 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document}α decay position.


Introduction
Experiments searching for rare events, such as neutrinoless double-beta (0νββ) decay [1], demand for a detailed background understanding in order to implement possible a e-mail: davide.chiesa@mib.infn.it (corresponding author) b e-mail: lorenzo.pagnanini@lngs.infn.it c e-mail: stefano.pozzi@mib.infn.it 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νββ 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 α-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 × 10 −3 counts/ (keV kg yr) in the region of interest around the 82 Se ββ decay Q value (Q ββ = 2997.9 ± 0.3 keV [13]), characterized by an average energy resolution of (20.05 ± 0.34) keV FWHM [14,15]. Such impressive low background rate was achieved by combining the α-particle identification (and rejection) with the analysis of time-delayed coincidences. In particular, we tagged potential 212 Bi α decays and we vetoed any event occurring within 7 half-lives of its daughter 208 Tl (T 1/2 = 3.05 min), that β decays with a high Q value (5 MeV). In this way, we reduced the background in the region of interest by a factor ∼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 α-decay sequences in 232 Th and 238 U chains to improve the CUPID-0 background model [18]. Indeed, the α-decay features, especially the short range of energy deposition, allow to study the contaminant localization. In this paper, we describe how we analyzed the α − α 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νββ decay, because the background index induced by such contaminations strongly depends on the their location.

Experimental setup
CUPID-0 is the first 10 kg-scale CUPID [2] demonstrator using enriched scintillating calorimeters to search for 0νββ 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 µm thick Vikuiti TM reflective foil to enhance light collection. To measure the light signal, 170 µ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 ∼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].

Data production
In this work, we analyze the spectrum of α-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 α-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 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 the α-particles relying on the light pulse shape parameter defined in Ref. [24], which allows to discriminate >99.9% of α events from β/γ ones at energies >2 MeV. We tag the events that simultaneously trigger more than one crystal within a ±20 ms time window, assigning a multiplicity label (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 (∼ 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.

Search for delayed coincidences
In CUPID-0, the α-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 (Δt w ), determines the probability of random coincidences: Since P random has to be kept 1 and Δ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 1/r only.
In CUPID-0, the detector time resolution is of the order of few ms and the rate of α events is at maximum 1.7 × 10 −4 Hz/crystal and, on average, 6.3 × 10 −5 Hz/crystal. Therefore, the most suitable α-decay sequences in 238  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 (M 1 ) not piled-up α-events at the Q value peak of the first decay in the sequence, within a ±1.5 σ 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 Δ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 (∼97%), while keeping random coincidences at a negligible level. We only require that a daughter is an α-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 pileup 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 Δ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 Δt between couples of parent-daughter events. The expected number of these random coincidences between parent candidates is given by: where N ch P is the number of candidate parent events detected by each channel ch, and r ch P 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 P P ) 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 M 1 α-events passing the pile-up rejection cut (M 1α ).
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 ∼1% of the M 1α spectrum. To determine the number of detected delayed coincidences (N C ), i.e. couples of timecorrelated parent-daughter events, we exploit the daughter energy signature and we count the number of events falling in a ±3σ 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 α-α piledup event, that provides an unambiguous signature to identify this sub-chain of 3 consecutive α decays. We conservatively discard the piled-up events spaced less than 80 ms in time 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 α-decays, the third occurring with a Δt > 80 ms because above this threshold we are able to precisely trace back the second pulse amplitude to a full-energy 216 Po decay deposition.
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 Δt between couples of parent-daughter events. Indeed, the Δ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 Δt distribution of random coincidences can be approximated as flat when Δt w 1/r . As shown in Fig. 3, the measured Δ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: where r ch D is the α-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

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 Fig. 4 Sketch of α −α 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 α emitted in the daughter decay calorimeters consists in analyzing the spectrum of M 2 events comprised of α-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 α-particles escaping from their surfaces, preventing the detection of α-recoil M 2 events. To overcome such limitation, we conceived an innovative method based on the analysis of α-α 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 D Q |P Q , is almost 1. Conversely, when contaminations are on crystal surfaces, the α-particles have a not negligible probability to escape the detector. Thus in this case, the conditional probability p D Q |P Q to detect a daughter event at the Q value is significantly < 1.
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: 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 (ε). In this system, the different value of p D Q |P Q exploited to disentangle bulk and surface contaminations (hereafter labeled as p C ) enters in the ε C terms, which can be expressed as ε C = ε P p C . By solving the system in Eq. (3), we finally obtain the formula to calculate r : The physical constraint to be respected in order to get positive results for r is: 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.

Evaluation of delayed coincidence probability
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 (ε b P and ε s P ) and the probabilities to detect delayed coincidences ( p b C and p s C ) 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 ±5 ms window, and we group into multiplets the events involving different crystals within ±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/λ , where λ is the so-called depth parameter. This parameter is usually chosen in the range from few nm up to tens of μ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 λ values equal to 10 nm and 10 μm, being much lower and on the same scale of the α-particle range, respectively. Even for very shallow surface contaminations, we get a significant fraction of events reconstructed at the α-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 μ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 λ = 10 nm to model surface contaminations and to study the ratio between bulk and surface activities.
In Table 2, we report the values of ε b P , ε s P , p b C , and p s C obtained from the MC simulations of 238 U and 232 Th decay chains. The ε P efficiencies are computed by taking into account that a ±1.5 σ range was used to select the candidate parents in the experimental data. As expected from a simple geometric reasoning about α 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 α-decay  decreases as the contaminants are simulated in a thinner surface layer due to α-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 (λ = 10 nm) sequence and we have to discard a fraction of piled-up events with Δt < 80 ms to be consistent with the experimental data processing.

Results and discussion
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 ∼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 σ 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 α-decay sequences for the 238 U sub-chain: 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).

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 αs and nuclear recoils. The α escape probability is significantly affected when the contamination depth is changed from λ = 10 nm to a value on the same scale of the α range. For example, by analyzing our data with λ = 10 μm, the activity ratios r would scale up by a factor ∼ 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 λ 10 nm. Since in the experimental spectrum there are no emerging peaks at the α energies of the isotopes analyzed in this work, we can exclude that shallower contaminations (λ 10 nm) can significantly affect our results.

Conclusions
In this work, we presented an innovative analysis technique to study the background sources in cryogenic calorimeters relying on the time-correlation of α-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 coinci-dences in CUPID-0 allowed to measure the half-life of 216 Po [38].