Modeling of GERDA Phase II data

The GERmanium Detector Array (Gerda) experiment at the Gran Sasso underground laboratory (LNGS) of INFN is searching for neutrinoless double-beta (0νββ) decay of 76Ge. The technological challenge of Gerda is to operate in a “background-free” regime in the region of interest (ROI) after analysis cuts for the full 100 kg·yr target exposure of the experiment. A careful modeling and decomposition of the full-range energy spectrum is essential to predict the shape and composition of events in the ROI around Qββ for the 0νββ search, to extract a precise measurement of the half-life of the double-beta decay mode with neutrinos (2νββ) and in order to identify the location of residual impurities. The latter will permit future experiments to build strategies in order to further lower the background and achieve even better sensitivities. In this article the background decomposition prior to analysis cuts is presented for Gerda Phase II. The background model fit yields a flat spectrum in the ROI with a background index (BI) of 16.04−0.85+0.78·10−3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {16.04}_{-0.85}^{+0.78}\cdotp {10}^{-3} $$\end{document} cts/(keV·kg·yr) for the enriched BEGe data set and 14.68−0.52+0.47·10−3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {14.68}_{-0.52}^{+0.47}\cdotp {10}^{-3} $$\end{document} cts/(keV·kg·yr) for the enriched coaxial data set. These values are similar to the one of Phase I despite a much larger number of detectors and hence radioactive hardware components.


Introduction
A large fraction of current experimental efforts are devoted to test the precision of the standard model of particle physics and investigate the presence of new phenomena. Many extensions of the standard model predict rare processes and in particular the existence of neutrinoless double-beta (0νββ) decay [1][2][3]. The observation of this lepton-number violating decay would shed light on the nature of neutrinos and could give a hint on the scale of neutrino masses. The GERmanium Detector Array (Gerda) experiment [4,5] is searching for the 0νββ decay of the candidate isotope 76 Ge at a Q-value of Q ββ = 2039.061(7) keV [6]. Gerda is operating 37 detectors made from material enriched in 76 Ge and a total mass of 35.6 kg bare in 64 m 3 of liquid Argon (LAr, purity 5.0). The experiment profits from the high shielding JHEP03(2020)139 power of the LAr and its scintillation properties. A hybrid instrumentation consisting of light guiding fibers read out by silicon photomultipliers (SiPMs) and 16 photomultipliers (PMTs) detect LAr scintillation light in order to veto events depositing energy in the cryogenic liquid [5]. The LAr cryostat itself is situated inside a tank filled with 590 m 3 of purified water shielding against external ionizing radiation and neutrons. Furthermore, it is instrumented with 66 PMTs to veto muons by the detection ofČerenkov light. Gerda is the first 0νββ decay experiment working in a "background-free" regime in the region of interest (ROI) after analysis cuts [7][8][9], where the ROI is Q ββ ± FWHM/2, and FWHM is defined as full width half maximum.
In the following, we present the spectral decomposition of Gerda Phase II data. The analysis is conducted prior the application of active background suppression techniques to data, i.e. the LAr veto [5] and pulse shape discrimination (PSD) taking advantage of particular detector signal shapes [10]. A new assay of the Gerda background is necessary due to substantial upgrade works finished in 2015 [5]. Most structural components close to the detectors have been exchanged using materials with improved radio-purity, the detector array has been enlarged and the LAr veto instrumentation has been deployed during the upgrade. Moreover, each detector string (enclosed in a copper mini-shroud during Phase I) has been encapsulated in a transparent nylon mini-shroud in order to limit the drift of 42 K ions in the detector vicinity and appropriately propagate the LAr scintillation light [11] (see section 2.3 for details). The introduction of these new setup components and materials changes the distribution and composition of radioactive impurities in the setup.
A precise knowledge of the spectral composition of the data is a key point for further analysis like accessing the half-life of the lepton number conserving mode of double-beta (2νββ) decay. Moreover, there are significant efforts towards reaching the tonne-scale of active isotope mass and the localization of remaining radioactive impurities inside the setup is the basis for the possible further reduction of background. This is essential for future endeavors in order to boost the current signal discovery and limit setting sensitivity by two orders of magnitude to the range of T 0ν 1/2 > 1 · 10 28 yr.

Data selection and prior knowledge
The data analyzed in the following were taken between December 2015 and April 2018. In this period the Gerda array consisted of 40 high-purity germanium (HPGe) detectors: 30 Broad Energy Germanium (BEGe) detectors [12,13] and 10 detectors with a semi-coaxial geometry three of which are made from germanium with a natural isotope composition. The enrichment fraction of the 30 enriched BEGe ( enr BEGe) detectors is 87.8 % while the respective fraction for the 7 enriched coaxial ( enr Coax) detectors is in the range of 85.5-88.3 % [5].
Detector geometries. The Gerda HPGe detectors are made of p-type germanium. p + and n + contacts are manufactured via boron implantation and lithium diffusion, respectively. The p + electrode is connected to a charge sensitive amplifier while the n + electrode is biased at typically 4 kV. A groove between the two contacts provides electrical insulation. The bias high-voltage creates an internal electrical field which is responsible for charge JHEP03(2020)139 collection. When biased at full-depletion voltage, the germanium detectors reach maximal (ε = 1) charge collection efficiency (CCE) in an internal active volume, surrounded by a transition layer (TL) with reduced CCE (0 < ε < 1) and low electric field. The TL is covered by a thin conductive layer in which all charges recombine and charge collection is entirely suppressed (ε = 0), therefore, called dead layer. We define the contact thickness as the depth at which the CCE reaches its maximal value. The Gerda detectors are of two distinct geometries. In the semi-coaxial layout the thin p + contact (0.5 − 1 µm) covers the entire bore hole; in the BEGe-type, instead, the same contact is a disk of 15 mm diameter (see figure 3 in reference [14]). The n + contact, about 1 mm thick, "wraps around" the detector. An exhaustive description of the Gerda detector geometries and properties can be found in previous publications [5,[12][13][14]. The detector arrangement in the 7 strings that constitute the Gerda array is graphically presented in figure 1a (and in the appendix in figure 8).
Data acquisition and treatment. All data are recorded using FADCs and are digitally processed off-line [5]. The linearity of the data acquisition system and off-line energy reconstruction was tested with a precision pulse generator over the whole dynamic range of the FADCs. Up to an energy of at least 6 MeV no major non-linearity and pulse shape deformation was observed. A signal above threshold in any of the germanium detectors triggers the data acquisition and the respective event is written to disk. An event is defined as the set of traces recorded in the 40 germanium detectors, 16 photomultipliers (PMT) and 15 silicon photomultiplier (SiPM) channels from the LAr veto and the signal from the WaterČerenkov muon veto. This hardware threshold is detector and run dependent and varies between 20 keV and 200 keV. The energies of all other detectors are reconstructed from the recorded traces and we apply a threshold of 40 keV on these. At this threshold the reconstruction efficiency is practically 100 %. Hence, we are sure to treat data and simulations in a consistent manner. In the following, we define the multiplicity of an event as the number of germanium detectors in which an energy of at least 40 keV is registered.
The energy deposition associated to each germanium detector signal is determined via a zero area cusp (ZAC) filter [16] which is optimized off-line for each detector and each calibration. Calibrations are usually taken with three 228 Th sources which are lowered into the LAr to the vicinity of the detector array in a 1-2 week cycle. An energy correction due to crosstalk between detector channels is performed for each event. The average crosstalk for all pairs of channels is about 0.05%. Details about the crosstalk correction can be found in reference [17]. Events in a window Q ββ ± 25 keV are excluded from the analysis until all selection cuts are finalized. The number of events and their energies in this window are only released once all analysis steps are defined.
Each event has to pass a number of quality cuts which are tailored to filter unphysical events [7]. Data taking periods in which stable operation cannot be guaranteed are excluded from analysis. The overall duty cycle in Gerda Phase II is 92.9 %. We keep 80.4 % of the recorded data as valid analysis data set, discarding for instance, periods of unstable calibration or earthquakes like in August 2016. Detectors with an unstable energy calibra-  [13]. Hence, events triggered by this detector are not considered in either data set. GD02D is omitted from the mass computation. Table 1. Properties of the data sets considered in this analysis. Further details about the Gerda detectors can be found in past publications [13,14].
tion are used only to determine the event multiplicity but do not enter any data set, e.g. an event that triggers three detectors one of which cannot be calibrated well is not considered a two-but a three-detector event. Also, two-detector events involving a detector which is not well calibrated are rejected. Events with a multiplicity higher than two are discarded by default and, likewise, events which trigger the muon veto are excluded.

Analysis data sets
Events of multiplicity one (M1) and multiplicity two (M2) from detectors with enriched isotope composition are accounted for in the construction of the analysis data sets. Events from the coaxial detectors with natural isotope composition, located in the central detector string, are not used in this analysis due to large uncertainties on their n + contact thickness and detection efficiency. The M1 events are split in two data sets based on the two enriched detector geometries which we call M1-enrBEGe and M1-enrCoax in the following. The M2 data form a third data set which is named M2-enrGe. The energy we associate to an M2 event is the sum of the energies reconstructed in the two detectors. The data sets, their exposure and respective detector mass are listed in table 1.

Monte Carlo simulations and probability density functions
The Probability Density Functions (PDFs) used to model contributions to the energy spectra are obtained from Monte Carlo simulations. The latter are performed using the MaGe simulation framework [18], based on Geant4 v10.4 [19][20][21]. MaGe contains a software implementation of the Gerda Phase II detectors as well as the assembly and all other surrounding hardware components. A visualization of this implementation is presented in figure 1. Detector intrinsic 2νββ decays of 76 Ge and background events originating from radioactive contaminations in and around the detector assembly are simulated. The energy spectrum of the two electrons emitted in the 2νββ decay was sampled according to the distribution given in reference [22] implemented in Decay0 [23]. The PDFs are obtained from the Monte Carlo simulations, taking into account the finite energy resolution and individual exposure acquired with each detector during the considered data taking periods. Special care is taken not to statistically bias the PDFs by assuring that each simulated decay is taken into account only once in the production of a PDF. For more details see appendix C. JHEP03(2020)139 (e) (f) Figure 1. Implementation of the Gerda array in MaGe, visualized using the Geant4 visualization drivers. From left to right: a) the Gerda detectors, b) the holder mounting, composed of silicon plates and copper bars c) the high-voltage and signal flexible flat cables plus the front-end electronics on top, d) the full array instrumentation, including the transparent nylon mini-shrouds, e) the full LAr veto system surrounding the array, including the fiber shroud (in green), the Tetratex Rcoated copper shrouds (above and below the fibers) and the two PMT arrays, f) the LAr veto system without the copper shrouds.

Background expectation
The event energy distribution of the three data sets is displayed in figure 2; the sum spectrum of M1-enrBEGe and M1-enrCoax in the top panel and M2-enrGe in the bottom panel. For the single-detector data, in the top panel, the following features are most noticeable: the β decay of 39 Ar dominates the spectrum up to 565 keV while between 600 and 1500 keV the most prominent component is the continuous spectrum of 2νββ decay of 76 Ge. Two γ lines at 1461 and 1525 keV can be attributed to 40 K and 42 K; further visible γ lines belonging to 85 Kr, 208 Tl, 214 Bi and 228 Ac are indicated in the figure. The highest energies displayed are dominated by a peak like structure emerging at 5.3 MeV with a pronounced low energy tail. This is a typical spectral feature of α particles and can, here, be attributed to 210 Po decay on the thin detector p + surfaces [14]. Events above the 210 Po peak belong to α decays emerging from the 226 Ra sub-chain on the detector p + surfaces. All these components contribute also to M2-enrGe except for 39 Ar, 2νββ and high energy α components. This is due to the short range of α (tens of µm) and β particles (typically smaller than 1.5 cm) in LAr and germanium with respect to the distance between detectors which is of the order of several cm. The structural components of the setup have been screened for their radio-purity before deployment. Two measurement methods were used depending on the screened isotope: γ ray spectroscopy (Ge-γ) with High Purity Germanium (in four underground laboratories, for details see reference [4]) and mass spectrometry with Inductively Coupled Plasma Mass Spectrometers (ICP-MS) [24]. Especially materials close to the detectors have been screened for radioactive contaminations originating from the 238 U and 232 Th decay chains, 40 K and 60 Co. For measured specific activities and upper limits see reference [5] section 5. All possible background sources taken into consideration in this analysis are described in detail below. The descriptions are accompanied by a selection of PDFs in figure 3 (see also appendix C). Their decay products consist of γ or β particles with an energy higher than 520 keV. Less energetic particles from the remaining constituents in the chain do not enter the energy window which is considered in the presented analysis. The α emitters from the decay chains contaminating the thin p + electrodes are described below.

JHEP03(2020)139
60 Co. A significant fraction of components in the Gerda setup is made of copper [5], which can be produced with high radio-purity but is potentially activated by cosmic rays producing the long-lived isotope 60 Co. The latter decays with a half-life of 5.2711(8) yr; from material screening it is also expected to be found in some of the detector high-voltage flexible flat cables. 40 K. This isotope is found in all screened materials. Construction materials were not optimized for ultra-low 40 K content because the Q-value of its decay is well below Q ββ and hence does not contribute to the background in the ROI. The 40 K decay spectrum exhibits a γ line at 1460.822(6) keV with an accumulated statistics on the order of 100 cts/detector. In figure 12 the expected counts per detector for 40 K simulated in different locations are shown. Using the ratio of events detected in different detectors, information about the spatial distribution of 40 K can be extracted. We use this spatial information to resolve degeneracies of 40 K in the energy spectra (for details see appendix A). 42 K. A cosmogenically produced isotope in LAr is 42 Ar (T 1/2 = 32.9(11) yr) which decays to 42 K. The distribution of 42 K inside the LAr is likely to be inhomogeneous due to drift of the ionized decay products induced by the electric field (generated by high-voltage cables and detectors) and convection. 42 K decays to 42 Ca via β decay with a half-life of 12.355(7) h and a Q-value of 3525.22 (18) keV, well above Q ββ . For the β particle to be detected the decay needs to happen within a distance of a few centimeters 1 to the detector surface. As the detectors are in direct contact with the LAr, the β component of 42 K potentially gives one of the most significant contributions to the background in the ROI. Therefore, we separate decays originating inside and outside the mini-shrouds in the following analysis. The full-range fit has little sensitivity to any potassium inhomogeneity outside the minishrouds. Based on detector-wise observations, however, a surplus of 42 K above the detector array in the vicinity of the front-end electronics is deduced (see appendix A). Outside the mini-shrouds we, hence, consider a homogeneous component and an additional distribution above the detector array. Inside the mini-shrouds the β spectrum becomes potentially important. Some scenarios are possible, the closer 42 K decays to the detector surface, namely to the n + and p + contacts, the more β particles enter the germanium. A fraction of events around Q ββ coming from 42 K is potentially due to γ particles with higher energy and sub-percent level branching ratio or simultaneous energy deposition of multiple γ particles. This γ component could become important for large quantities of 42 K not located directly on the detector surfaces with the β particle being absorbed in the LAr. As for 40 K also the γ line at 1525 keV of 42 K contains valuable information about the spatial decay distribution of this isotope. In contrast to 40 K no additional information, e.g. from radiopurity screening measurements, is available. For more detailed information about 40 K and 42 K see appendix A.

JHEP03(2020)139
α emitters. The lithium-diffused n + detector surfaces act as a barrier for α particles. The latter can only penetrate the very thin boron-implanted p + -contact or the contact separating groove. α particles have to be emitted directly at the surface or from a thin adjacent layer of LAr. Since α particles have to cross the ∼ 0.5 µm thick p + dead layer and therefore only part of their initial energy is deposited in the active volume, this background component leads to peaks with characteristic low-energy tails in the HPGe energy spectra (see figure 3e). Some α events, presumably originating from the detector groove, are reconstructed with degraded energy and lead to an additional, continuous spectral component. We find mainly 210 Po but also traces of isotopes from the 226 Ra decay chain.
Detector bulk impurities. Cosmogenically produced long-lived isotopes can also be found in germanium [25][26][27]. In particular, 68 Ge and 60 Co can occur as detector intrinsic impurities with half-lives of 270.93(13) d and 5.2711 (8) yr. The BEGe detectors were kept underground during major parts of the fabrication and characterization operations. Periods when these detectors were above ground have been tracked in a database [12]. Thus, for the well-monitored BEGe detectors we expect impurities of 5 nuclei/kg of 68 Ge and 21 nuclei/kg of 60 Co as of September 2014 [12]. Extrapolating the expected impurities to the whole Phase II data taking period we expect on average 0.03 cts/day from 68 Ge and 0.1 cts/day due to 60 Co. From background modeling in Phase I [14] the contribution for the coaxial detectors formerly used in the Heidelberg-Moscow (HdM) [28] and Igex [29] experiments is expected to be even smaller due to their long storage underground. Simulating the expected detector bulk impurities we find background contributions around Q ββ of less than 10 −4 cts/(keV·kg·yr) in both cases. Hence, we conclude that 68 Ge as well as 60 Co can be neglected in the following analysis. Potential bulk contaminations with 238 U and 232 Th were studied in reference [30]. Only upper limits were found, establishing germanium crystals as material of outstanding radio-purity. Hence, we only consider the decay of 76 Ge via 2νββ as detector intrinsic background component while all other intrinsic impurities are considered to be negligible.
Other sources. As discussed in reference [14], prompt cosmic muon induced background events are efficiently vetoed by the identification ofČerenkov light emitted by muons when they pass the water tank. The expected BIs, due to the direct muon and neutron fluxes at the LNGS underground laboratory, have been estimated to be of the order 3 · 10 −5 cts/(keV·kg·yr) [31] and 10 −5 cts/(keV·kg·yr) [27] in earlier works, respectively. Background contributions coming from delayed decays of 77 Ge and 77m Ge, also induced by cosmic muons, are estimated to be 0.21 ± 0.01 nuclei/(kg·yr) [32] corresponding to a BI prior to the active background suppression techniques of about 10 −5 cts/(keV·kg·yr). Also, the water tank and LAr cryostat contaminations are expected to contribute to the Gerda BI with less than 10 −4 cts/(keV·kg·yr) [4,33]. All above mentioned contributions are considered negligible in this work. Other potential sources of background from interactions of 76 Ge [15,27] and 206 Pb [34] with neutrons and 56 Co for which no evidence was found are not taken into consideration. The cosmogenically produced isotope 39 Ar and the anthropogenic isotope 85 Kr [35], which are dissolved in LAr, emit particles which are dominantly less energetic than the energy window which is considered in the presented analysis.
(c) 40 K contamination close to the detector array (on the mini-shrouds), at a higher radial distance (on the fiber shroud) and higher vertical distance (on the copper shrouds).

JHEP03(2020)139 3 Statistical analysis
The multivariate statistical analysis, which is used to model and disentangle the background in its components, runs on the three binned data sets M1-enrBEGe, M1-enrCoax and M2-enrGe. It is based on the reconstructed energy with the zero area cusp (ZAC) filter algorithm which is close to optimal and provides an excellent low-frequency rejection [16]. The single-detector data sets M1-enrBEGe and M1-enrCoax contain the reconstructed ZAC energy of all M1 events whereas for the two-detector events the sum of the two reconstructed energies is put in the M2-enrGe data set. Moreover, the count rate per detector is used for the two potassium γ lines. The spatial event distribution is a collection of the number of events per detector for M1 events and expressed in a matrix of pairs of detectors for all M2 events.
Assuming that the number of events in each bin follows the Poisson probability distribution Pois(n; ν), where ν is the expected mean and n is the experimentally measured number of counts, the likelihood function for a binned data set reads N bins i=1 Pois(n i ; ν i ).
is the expected number of events in the i-th bin, calculated as the sum of the contributions from each background component k; ν i (λ 1 , . . . , λ Ncom ) is a function of the parameters of interests λ j (isotope activities, 2νββ half-life, etc.). The complete likelihood function adopted for the present analysis combines the three data sets M1-enrBEGe, M1-enrCoax and M2-enrGe: The statistical inference is made within a Bayesian framework. Hence, to obtain posterior probabilities for the free parameters of interest λ j , the likelihood defined in equation (3.1) is multiplied according to the Bayes theorem by a factor modeling the prior knowledge of each background component as presented in section 2.3. The computation is performed using a Markov Chain Monte Carlo (MCMC) and is implemented using the BAT software suite [36,37]. Posterior probability distributions of any observable that is not a free parameter of the likelihood function, like background index estimates, are obtained by sampling the desired parameter from the MCMC. A p-value estimate is provided as a goodness-of-fit measure by adopting the algorithm suggested in reference [38] for Poissondistributed data. It has to be kept in mind that this p-value estimate, however, is not as well suited for model comparison as is for instance a Bayes factor; e.g. the number of free parameters is not taken into account while a Bayes factor always penalizes models that add extra complexity without being required by the data.

Analysis window and binning
The fit range and data bins are chosen such as to exploit as much information from spectral features as possible brought by data without introducing undesired bias. The chosen fit range in energy space for the single-detector data sets (M1-enrBEGe and M1-enrCoax) starts from just above the end-point of the 39 Ar β − spectrum at 565 keV and ends just above the 210 Po peak at 5260 keV, where the event rate drops to almost zero values. For the JHEP03(2020)139 two-detector events (M2-enrGe data set) the fit range starts at 520 keV and extends up to 3500 keV. Possible additional components outside of this range (e.g. 39 Ar) do neither add information to the background decomposition in the ROI around Q ββ nor to the analysis of 2νββ decay. Furthermore, at energies lower than ∼100 keV the shape of the PDFs is dominated by uncertainties on the detector transition layer model, which describes the charge-carrier collection at the interface between the n + contact and the detector active volume. The exact nature of this transition region is different for each detector and prone to systematic uncertainties [39].
With an energy resolution which is typically 3-4 keV at Q ββ (FWHM) [8,9] and better at lower energies, a fixed bin size of 1 keV was chosen for all data sets. The only exceptions are the two γ lines from 40 K and 42 K each of which is combined in a single bin from 1455 keV to 1465 keV and from 1520 keV to 1530 keV, respectively. This is done in order to suppress any systematic uncertainties of the energy calibration and resolution model that affect the position and shape of the γ lines [9].

Likelihood factorization
A feature of the selected data is that the likelihood in equation (3.1) can be factorized in uncorrelated parts which can be studied individually and in detail. In the following we shortly outline the parts of the data which were studied in depth based on the approach of factorizing the likelihood into uncorrelated parts. Finally, the results of these analyses are incorporated into a full-range fit. This procedure is equivalent to a simultaneous analysis of all data but increases the input knowledge for the fit and breaks down the computational complexity in smaller steps.

JHEP03(2020)139
α events background analysis. The single-detector energy spectra above 3.5 MeV (the Q-value of 42 K β decay) are strongly dominated by α events. They are not present in twodetector data due to the short range of α particles in LAr and germanium. Also, this component is not correlated to other backgrounds considered here because it peaks at energies well above the highest γ emission energies and β decay Q-values. A careful study was carried out considering various p + contact thickness and event rates to reproduce the 210 Po peak. In order to reproduce α events with degraded energy an empirical model is fit to the data. A linear function with free slope and offset and a cut-off below the maximum of the 210 Po peak fits the data well. The agreement of the α background model with the data is demonstrated in appendix B and figure 9 therein. Information from the detailed analysis of the high-energy α region is incorporated in the full-range fit using a combined PDF that summarizes the 210 Po peak plus the 226 Ra decay chain and a linear floating component for degraded α events.

Prior distributions
The following criteria are adopted to convert the prior information described in section 2.3 into prior probability distributions on the parameters of interest: 3 if a measured value with uncertainty is available for a background contamination then a Gaussian distribution with a corresponding centroid and a 1σ width is adopted. In presence of a 90 % C.L. upper limit, instead, an exponential prior distribution is constructed with 90 % of its area covering parameter values from 0 up to the given 90 % C.L. upper limit. A uniform prior distribution is assigned to components for which no measured value or upper limit is available. Ranges for uniform priors are initially taken very wide, in order to span a large portion of the allowed parameter space, then optimized to contain at least 99 % of the posterior distribution. As mentioned before, in addition to the information from screening measurements, prior distributions for 40 K and 42 K are constructed considering the posterior inference from their spatial distribution. 4 Moreover, as 214 Bi is part of the 226 Ra decay chain, we constrain a 214 Bi component on the p + contact by a Gaussian prior extracted from the obtained 226 Ra activity based on the energy estimator in the high-energy α region.

Results
As described in section 3.2 the α event background and potassium γ lines are studied individually and the results are incorporated in the full-range fit as prior distributions. The latter combines a simultaneous fit of the M1 and the M2 data sets. For the final combination of parameters, outlined in this section, components with a posterior distribution peaked at zero were eliminated from the fit. The stability of the results with respect to the bin size and prior distributions was verified. Changing the prior distribution for fit parameters for which no screening measurement is available from a flat to an exponential one does not 3 In Bayesian analysis the prior probability distribution describes all knowledge about an unobserved quantity of ultimate interest before taking the data into account. 4 The Bayesian posterior distribution is the conditional probability distribution of the unobserved quantities of ultimate interest, given the observed data.

JHEP03(2020)139
significantly impact the final posterior distributions. The compatibility of the final model, which includes 34 free fit parameters, with data is supported by a p-value of ∼ 0.3.
The estimated activities of individual components and other parameters of interest are listed in table 2. In particular, for each component we report the global and the marginalized mode of the posterior parameter distribution, along with its smallest 68 % C.I. The global mode corresponds to the global best fit value while the marginalized mode is the most probable parameter value when integrating over all other parameters. The original type of prior distribution is marked with [f] for flat, [g] for Gaussian and [e] for exponential; the latter two are used if screening measurements are available. Subsequently, for all 40 K and 42 K components, the prior distribution is imported from the potassium tracking analysis and for 214 Pb and 214 Bi on the p + contact from the reconstructed 226 Ra content from the α events background analysis.
The spectral decomposition of all data sets is shown in figure 4. For each data set the residual distribution as a multiple of the expected 1σ fluctuation in each bin is displayed. We find for the M1-enrBEGe data set 66.4 %, 94.5 % and 99.6 % of points in the 1σ-, 2σand 3σ-bands, for the M1-enrCoax data set 66.0 %, 94.7 % and 99.8 % and for the M2-enrGe data set 70.0 %, 96.1 % and 99.7 %, respectively. Thus, in all three cases the residuals are normally distributed. No outliers with residuals larger than 3σ are found in a ±50 keV window around Q ββ and the bins exceeding 3σ do not correspond to any noted γ line.
The 42 K distribution is optimized to best fit the data. In order to disentangle the 42 K γ and β components, the volume inside and outside of the mini-shrouds is separated in the PDF construction. Inside the mini-shrouds a homogeneous distribution is compatible with the data as well as 42 K attached to the detectors contact surfaces. In the fit model given here, a possible scenario is chosen where all 42 K is located on the n + surfaces. However, we note that 42 K on the p + appears to partly substitute the energy-degraded α component in the M1-enrCoax data set if introduced in the fit and predicts a higher total BI around Q ββ . The extracted 42 K activity on the enr Coax p + contact in this case is 22 ± 4 µBq corresponding to a contribution to the BI around Q ββ of (7 ± 1) · 10 −3 cts/(keV·kg·yr). For the M1-enrBEGe data set the posterior distribution of a possible 42 K component on the p + contact is compatible with zero. Outside the mini-shrouds an inhomogeneous distribution of the 42 K decays explains the observations better. Detectors which are located at higher positions in the strings show an excess of events in the 42 K 1525 keV γ line which is compatible with a surplus of 42 K located right above the detector array (see appendix A). The full-range fit model contains a homogeneous 42 K distribution outside the mini-shrouds which is reconstructed with a specific activity of 186 ± 39 µBq/kg plus an additional distribution in the vicinity of the cables above the array.
A large fraction of the contamination with 40 K in the setup cannot be accounted for by the screened hardware listed in table 2. We thus add a close (∼ 1 cm) and a far (∼ 50 cm) 40 K component with respect to the detector array which are in fact replica of the PDFs for the mini-shrouds and the Tetratex R -coated copper shrouds. These additional components absorb the excess indicated by the fit, the largest part of the reconstructed events in the spectra is attributed to impurities close to the array.   Note that for bins with low expected statistics due to the discrete nature of the measured spectrum not all colored bands are meaningful [40].
The 40 K and 42 K distributions can be further split into smaller volumes and studied as an extension of the potassium tracking analysis (as described in section 3.2) projected in detector space. The additional 40 K component close to the array and the 42 K component above the array are split into 7 sub-components on a string-by-string basis. The potas-JHEP03(2020)139 sium concentration is in general found to be asymmetric among the detector strings. In particular, a more prominent 42 K concentration is found above the central string. This is consistent with the electrostatic drift of 42 K ions induced by the electric field in the LAr which is generated by the unshielded high-voltage flat cables biased with about 4 kV. The 40 K and 42 K spatial analysis fitting the potassium γ lines projected in detector space is presented in full detail in appendix A.
The α distribution is adjusted to best fit the data. The 210 Po peak at 5.2 MeV is found to be best described by a mixture of PDFs obtained assuming different p + contact thicknesses confirming results of the Phase I background analysis [14]. The empirical linear model which is used to describe α events with degraded energy (see section 3.2), extends down to Q ββ and below. For the M1-enrBEGe data set α events are efficiently isolated using pulse shape discrimination (PSD) techniques. The compatibility of the degraded-energy α component with α events identified by PSD was checked and is found consistent. All details about the α events analysis can be found in appendix B.
Smaller contributions to the background model in the full energy range are attributed to 214 Pb and 214 Bi from the 238 U decay chain, 228 Ac, 212 Bi and 208 Tl from the 232 Th decay chains and 60 Co. With a total contribution in the fit range of 10 −3 cts/keV for both the M1-enrBEGe and M1-enrCoax data set 234m Pa gives negligible contribution to the spectra and is therefore dropped from the full-range fit model. The central values preferred in the full-range fit are driven by screening measurements and the spectral contributions are all fully accounted for by the listed hardware components. The only exception is 214 Pb and 214 Bi where a minor contribution is added on the p + contact expected from the observation of α events belonging to the 226 Ra decay chain.
Most counts in the fit range are attributed to the 2νββ decay of 76 Ge; in fact its continuous distribution dominates the spectrum up to almost 1.9 MeV. Here, we base the 2νββ half-life estimate on the M1-enrBEGe data set only. An additional parameter, δ 2ν , parameterizes the observed discrepancy to the value solely derived from the M1-enrCoax data set. The value of δ 2ν extracted from the fit amounts to a surplus of 5 % of 2νββ counts observed in M1-enrCoax. It mainly quantifies the systematic biases between the active volume determination methods of the two detector types. The enr BEGe detectors active volume measurements are affected by a smaller systematic uncertainty than the enr Coax detectors [13,14]. Hence, the extracted 2νββ half-life, based on the M1-enrBEGe data set and given here only with statistical uncertainties, amounts to T 2ν 1/2 = (2.03 ± 0.02) · 10 21 yr. A detailed discussion follows in section 5.
The background model describes the individual contributions to the total BI around Q ββ prior active background suppression (see figure 5). The BI is defined as the number of counts over exposure and energy in the energy window from 1930 keV to 2190 keV excluding the region around Q ββ (Q ββ ± 5 keV) and the intervals 2104 ± 5 keV and 2119 ± 5 keV, which correspond to known γ lines from 208 Tl and 214 Bi. The values for each background contribution are given in table 2. The dominating background contribution around Q ββ in the M1-enrBEGe data set come from 42 K. Isotopes from the 232 Th decay chain, α particles mainly with degraded energy and isotopes from the 238 U decay chain contribute about equally. The estimated total BIs extracted from the marginalized posterior distributions JHEP03(2020)139

Discussion
In general, impurities close to the detector array contribute most to the background, far components give minor contributions. The posterior distribution and the screening measurements are in very good agreement and the spectral content of each source of background can be accounted for by the screened hardware components. Only in the case of 40 K a large part of the observed activity cannot be explained by the screened hardware and is fit with the additionally introduced components far and close to the detector array. The 42 K and α event distributions cannot be constrained by screening measurements and are adjusted to best fit the data. The presented background model is not unambiguous in all components. As shown in figure 3 several PDFs of the same source of background located in different structural components are very similar and thus prone to correlation. Most of them have been resolved by introducing prior distributions based on the screening measurements. However, a few anti-correlations persist which are listed in table 3.

JHEP03(2020)139
For what concerns 42 K in the LAr volume outside the mini-shrouds and thus more distant from the detector array, the adopted distribution is purely empirical. Our prior knowledge is limited by the fact that the 42 K ions undergo drift due to the electrical fields surrounding the detectors and high-voltage cables. Also, due to thermal gradients they can be displaced by convection. Hence, their distribution inside the Gerda LAr is prone to systematic uncertainties. The presence of unshielded high-voltage cables above the detector array can explain the excess of 42 K found in this region. From the perspective of the full-range fit a more sophisticated modeling does not significantly modify the 42 K PDFs and hence the fit results. A potentially asymmetric 42 K distribution is, thus, not further followed in the main analysis. Nevertheless, some considerations can be found in appendix A. An explanation for 42 K on the p + contact being rejected for the M1-enrBEGe data set but potentially present in the M1-enrCoax data can be the specific bore-hole geometry of the semi-coaxial detectors. 42 K produced inside the hole cannot easily escape and is trapped close to the p + contact.
For each source of background the contribution to the BI at Q ββ prior to active background reduction is listed in table 2. The statistical uncertainties on the single contributions to the BI are generally of the order of 10 % or lower, with the exception of 42 K and energydegraded α events, for which the uncertainty is roughly doubled. The two contributions are affected by a higher uncertainty because they are not bound by screening measurements.
The background event distribution in the 0νββ analysis window can be well approximated with a constant function (see figure 5). With this assumption, the BIs extracted from data are 16.4 +1.7 −1.6 · 10 −3 cts/(keV·kg·yr) for M1-enrBEGe and 15.4 +1.8 −1.6 · 10 −3 cts/(keV·kg·yr) for the M1-enrCoax data set. These values agree well with the background model description presented in section 4. The BIs prior to further analysis cuts and before the upgrade of the Gerda experiment to Phase II can be found in reference [41]. For the M1-enrCoax data set the BI prior to the upgrade of (18 ± 2) · 10 −3 cts/(keV·kg·yr) is very consistent with the values presented here. The BI of the M1-enrBEGe data set instead is substantially improved from a Phase I value of 42 +10 −8 · 10 −3 cts/(keV·kg·yr) to a value which is at least 2.5× smaller in Phase II despite a significant increase of inactive hardware mass. 5 Contributions to the BI from all isotopes have been improved with respect to Phase I with the exception of background introduced by α surface events. The most drastic improvement is notable for 42 K for which the BI contribution for the enr BEGe detectors appears four times smaller than before the upgrade to Phase II.
As mentioned in section 4, the extracted 2νββ half-life estimate is based on the M1-enrBEGe data set only. The additional parameter δ 2ν mainly quantifies the systematic biases between the active volume determination methods of the two detector types. The full charge collection depth (FCCD), which determines the active volume of a detector, was studied extensively in a detector characterization campaign for the enr BEGe detectors [12,13]. The estimate of the FCCD used in this analysis is based on measurements using an 241 Am source with characteristic γ lines at 60 keV, 99 keV and 103 keV.

JHEP03(2020)139
However, the FCCD was also measured using a 60 Co source with characteristic γ energies of 1173 keV and 1332 keV. The latter FCCD Co is systematically higher (about 3 %) with respect to the FCCD Am . The discrepancy could be explained by an energy dependence of the initial charge-carrier cloud size inside the detector but the actual impact on the active volume is still under investigation. For the enr Coax detectors only FCCD values determined with a 60 Co source are available. Considering the systematic uncertainties affecting the determined active 76 Ge exposures of the M1-enrBEGe and M1-enrCoax data sets (1.8 % and 5 % respectively, see table 1) δ 2ν is compatible with zero within 1σ. 6 Various systematic effects have to be considered when estimating the uncertainty on the 2νββ half-life T 2ν 1/2 . Due to the fact that the aim of the paper is not a precise 2νββ half-life measurement, for most of them only a conservative evaluation is provided. Several systematic uncertainties arise from the Monte Carlo simulation framework. Uncertainties due to the Geant4 model of particle interactions and propagation were estimated to be of the order of 2 % in previous publications [42,43]. Approximations in the implementation of the Gerda setup are conservatively estimated within a 1-2 % uncertainty range. This accounts for possible spectral shape modifications due to inaccurate charge collection model between the n + contact layer and the active detector volume. Uncertainties induced by the theoretical model of 2νββ decays implemented in Decay0, as well as data acquisition and selection methods are considered negligible. A 1.8 % contribution accounts for uncertainties in the enrichment and active mass fraction determination (see active 76 Ge exposure in table 1). All the systematic effects considered above sum up to a total systematic uncertainty on T 2ν 1/2 of 3-4 %. In total this leads to T 2ν 1/2 = (2.03 ± 0.09) · 10 21 yr compatible with earlier results [42,43].

Conclusions
We presented the background decomposition of Gerda Phase II data before the application of active background suppression techniques using a multivariate Bayesian fit approach based on single-and two-detector data in energy and detector space. The model is able to well describe the data and the results are compatible with the expectations from material screening measurements. The only exception is 40 K for which a higher contamination is found, dominantly in hardware components close to the detector array. This indicates contaminations introduced during production and mounting procedures different from the screened reference samples; in fact a few parts underwent further processing after material screening. Analyzing the count rates in the 40 K and 42 K high-statistics γ lines on a detectorby-detector basis we find indications for asymmetries in the spatial distribution of the two potassium isotopes. Furthermore, the background indices at Q ββ prior active background suppression techniques are given by and are in very good agreement with the assumption of a flat background distribution in this region. In terms of the BI the upgrade to Gerda Phase II proves extremely successful. Despite major hardware changes and higher inactive mass close to the detectors, the BI before applying active background reduction remains unchanged for the enr Coax detectors and is improved by a factor of three for the enr BEGe detectors. A careful background model is essential in order to separate the two-neutrino doublebeta decay events from the other background components. We expect to substantially improve the precision of the T 2ν 1/2 measurement after applying the LAr veto cut. In this manner, the signal to background ratio in the 2νββ energy region is improved by about an order of magnitude [7,8]. Furthermore, this allows precision studies of the shape of the 2νββ spectrum and hence to test physics models beyond the Standard Model such as 0νββ decay with Majoron emission and Lorentz symmetry violation effects [43,44].
The localization of impurities makes the exchange of particularly contaminated components possible in upgrade works and thus the background can be potentially lowered even further. Moreover, it is important to learn what are the most important sources of background in order to improve handling and cleaning procedures as well as material selection. For future experiments like the Large Enriched Germanium Experiment for Neutrinoless ββ Decay (Legend) [45], which aims to cover the parameter space of inverted neutrino mass hierarchy, background reduction is the most crucial step in achieving the necessary sensitivity. The goal is to achieve a background index one order of magnitude lower than Gerda Phase II.

JHEP03(2020)139
modes to excited states of 2νββ decay in which the shape of the 2νββ decay spectrum is a unique feature and thus need to be well understood.
Initial observations in Phase II have shown that the 40 K and 42 K full-energy line intensities have increased by a factor of 4 and 2, respectively, in the single-detector data compared to Phase I [46]. The 42 K increase in activity can be attributed to the exchange of the mini-shrouds material from copper to nylon 7 during the Phase II upgrade: the electric field generated by the detectors bias high voltage is not screened by the conductive material anymore. The 42 K ions can be attracted from a larger LAr volume into the vicinity of the detectors. Moreover, the unshielded high-voltage cables could be an explanation for the higher rate of 42 K events seen in the uppermost detectors in the Gerda array. The higher 40 K event rate, on the other hand, is possibly attributable to the glue used for the nylon mini-shrouds and other new materials introduced with the LAr veto system. The exact amount, location and radio-purity of the glue is not precisely known. All changes to the setup that have been made during the upgrade to Phase II are described and motivated in exhaustive detail in reference [5].
In the following sections we focus on the characteristics of the events constituting the two potassium lines. In order to extract information about the spatial distribution of 40 K and 42 K contamination around the Gerda array, a treatment on a detector-by-detector basis is advantageous. The two γ lines contain enough statistics for such an analysis to be meaningful and constitute samples with a high signal to background ratio.

A.1 Data
Two windows around the potassium γ lines are projected in detector index space, such that, for single-detector data, each data point n i represents the total counts in detector i in the respective energy window. For two-detector data the detector space is two-dimensional, and each data point n ij represents the number of events for which energy is deposited in detector i and detector j.
The events in the potassium lines (denoted with K40 and K42 in the following) are selected in a ±3σ energy interval around the respective line, rounded up to an integer number of keV to match the specific energy windows in the energy distributions with 1 keV binning. σ is the energy resolution in the respective energy window. Additionally, three side-bands (SB1, SB2 and SB3 in the following) are used to estimate the continuum below and above the γ lines. Considering the further subdivision in single-(M1-) and two-detector (M2-) data, this leads to the definition of 5 × 2 energy regions, summarized in table 4. A visual representation of the selected windows can be found in figure 6. We use the PDFs respective to 214 Bi on the flat cables and detector intrinsic 2νββ decays to estimate the background. Other components are expected to contribute less in the respective energy windows.  Table 4. Energy ranges and corresponding number of events for the potassium tracking analysis (visualized in figure 6). Note that the windows for two-detector data are larger as the two singledetector energy resolutions are folded in the summed energy spectrum.

A.2 Analysis
The statistical approach of factorizing the likelihood is described in section 3.2. The part of the likelihood we are analyzing here runs simultaneously on the 5 × 2 energy ranges presented above. Following the naming convention introduced in section 3 it reads: where the index i runs over the bins (i.e. detectors) and the index d over the 5 energy windows considered, namely the three side-bands SB1, SB2, SB3 and the two line-bands K40 and K42. The M2-data sets are two-dimensional in detector space and run over the two indices j and k. Gaussian prior probability distributions for the 40 K activity are built from radio-purity screening measurements (see reference [5] section 5). For 42 K, for which no screening information is available, uniform priors are adopted, with the exception of the two 42 K components located on the n + contact surface of enr BEGe and enr Coax detectors. 42 K can be attracted to the n + surface by the electrical field created by the high voltage potential applied to the detectors. Both components are expected to be correlated by the volume ratio of the mini-shrouds (3:2 enr BEGe to enr Coax) the 42 K ions are attracted from. The volume ratio estimate is extracted from the geometric implementation in MaGe. We JHEP03(2020)139 assume an uncertainty of 0.1 mBq on either activity allowing for a change of their ratio. The correlation is included in the fit via a two-dimensional prior.
The analysis flow starts with a construction of a first, preliminary model, which consists only of background contributions that are expected from screening measurements of 40 K and known properties of 42 K. The resulting model, however, gives a non-satisfactory description of data and the posterior distributions for the 40 K components are significantly shifted to higher values with respect to the prior distributions, indicating a surplus of 40 K.
To find a better agreement with physics data while keeping the model as simple as possible, additional components using uniform priors are included one at a time in the fitting procedure, and the Bayes factor is calculated between the extended and the preliminary model. The model is iteratively updated by adding the component that results in the highest Bayes factor until no Bayes factor is larger than 10.
In a first iteration a replica of the PDF of 40 K in the mini-shrouds is added obtaining a Bayes factor 10. 40 K in the Tetratex R -coated copper shrouds is added in a second iteration with a Bayes factor of 11. For 42 K the only additional component that results in a Bayes factor greater than 1 is 42 K on the n + detector contacts. Although the fit shows only a slight preference (Bayes factor of 2) the component is added to the model because of its importance in the full-range fit, where the energy region above the 1525 keV γ line is also considered.
The results of the base model are shown in table 5 and a graphic representation showing the counts per detector in both potassium γ lines in M1-and M2-data can be found in figure 7. The analysis yields a p-value of ∼ 0.07, indicating an acceptable description of the data. To further improve the model rotationally asymmetric fit components are needed. The base model is accurate enough to be used as input for the full-range fit, which is insensitive to any rotational inhomogeneity of the location of background sources, as spectra from different detectors are merged into a single data set.
The two components 40 K close to the array and 42 K in LAr -above the array are split into 7 sub-components on a string-by-string basis (for the respective PDFs see appendix C). Furthermore, we consider a 40 K contamination on top of the central mini-shroud.
The results of this extended analysis are listed in table 6. An elevated 42 K concentration is found above the central string while a lower concentration is observed above the adjacent strings S1 and S6 (string numbers follow the nomenclature used in figure 8). Due to the large number of components the fit yields a high anti-correlation between the 42 K concentration above the outer strings and S7. This results in a high uncertainty on the latter fit parameter.
The screening measurements do not account for all observed 40 K. In general ICP-MS screening of the mini-shrouds with respect to 40 K is difficult and yielded only a lower limit. Different measurements seem to indicate different contamination levels of different mini-shrouds. Samples of glued nylon yielded the highest potassium contamination. As the gluing of the nylon mini-shrouds is done manually during installation the amount of glue and its exact location is hard to control. Hence, an asymmetric distribution is expected. The 40 K content of other close components like holders and cables might also be asymmetric. The asymmetric 40 K contamination is confirmed by the extended potassium JHEP03(2020)139 Figure 7. Decomposition of the energy windows corresponding to the two potassium lines in detector space: single-detector data (top) one-dimensional representation of two-detector data (bottom). Some components are merged for visualization purposes: in the K40 plots combined components are shown for 42 K and 214 Bi, while 40 K sources are grouped in close (flat cables, holders, mini-shrouds) and far (fibers, SiPMs, copper shrouds, front-end electronics) locations from the detector array. To visualize the two-detector data the sum of the projections on the two domain axes (index i and index j) is shown.

B α events background analysis
Above an energy of 3.5 MeV almost all registered events are due to α emitting isotopes. The respective part of the full likelihood can be approximately factorized and studied separately. α particles have a very short range in LAr as well as in germanium (continuous slowing down approximation, CSDA, range of 50 µm and 20 µm, respectively [47]) and are able to reach a detector's active volume only through the very thin (of the order of 500 nm) p + contact surface. Therefore, the α emitter contamination is detector-specific and depends only on the p + surface contaminations. Therefore, we analyze the enr BEGe and enr Coax detector data separately in energy space; the projection in detector space bares no correlation between detectors and hence contains no further useful information. The number of events in a single detector is not sufficient to further split the data on a detector-by-detector basis. The two data sets are uncorrelated and the statistical analysis can be carried out for each single-detector data set separately. In the two-detector data the α component is not observed due to the short range of these particles. All contaminations found are constituents of the 238 U decay chain. The main surface contamination observed is 210 Po which occurs either as an incident contamination and decays in time with a half-life of 138.3763(17) days [48] or is fed by a contamination with 210 Pb with a stable rate in time. The spectral form is identical for both cases and can only be disentangled by analyzing the α rate in time (see section B.1).
Above the 210 Po peak very few events are observed. In the M1-enrBEGe data set we find only four events with an energy larger than 5.3 MeV, while in the M1-enrCoax data set 22 such events are observed, 14 of which in a single detector ANG2 (see table 7). These events are due to α decays from 222 Rn and subsequent isotopes on the p + detector surfaces.
ANG2 also shows a higher 226 Ra (mother nucleus of 222 Rn) contamination which suggests dominantly a surface contamination with 226 Ra rather than 222 Rn dissolved in the LAr. In the latter case the decay chain would be broken as only the gaseous 222 Rn emanates from the construction materials into the LAr. The number of counts is too low to distinguish the spectral shape above 5.3 MeV and disentangle a surface contamination with 226 Ra from 222 Rn dissolved in LAr. A comparison between the counts observed above 5.3 MeV and the 214 Bi 609 keV γ line suggests that α events due to a dissolved 222 Rn contamination would not produce observable counts in said energy region. Assuming that all 214 Bi observed comes from dissolved 222 Rn leads, in fact, to a specific activity smaller than 10 µBq/kg. Hence, in the following, we will only consider a p + surface contamination with 226   Due to the very short range of α particles the energy spectrum of α decays exhibits a line with a pronounced low-energy tail. The tail is formed when the decay occurs under an incident angle with respect to the contact and the α particle loses part of its energy before reaching the detectors active volume. The maximum is shifted with respect to the full emission energy which is due to energy loss inside the electrode and depends on its minimal thickness. The detectors have slightly different contact thicknesses, also, the p + contact of a single detector may intrinsically be inhomogeneous. Therefore, we model the 210 Po peak with a mixture of PDFs obtained from simulations with different contact thicknesses. Due to the low number of counts observed in the 226 Ra chain it is sufficient to model this component with only one PDF. Furthermore, the isotope contamination is assumed to halve at each decay step. A reduction effect of the subsequent α decays in the 222 Rn chain had been observed in Phase I and attributed to possible recoil off the surface into the LAr [14]. We adopt this explanation in our model although we note that the number of events observed with an energy >5.3 MeV is not sufficient to confirm the previously observed reduction effect. Further details about the construction of the PDFs are given in appendix C.
Dedicated measurements [49] have shown that events originating in the contact separating groove are partly reconstructed with degraded energy. A simulation-based model of these energy-degraded events is not available yet. We approximate this component with an empirical linear distribution truncated below the maximum of the 210 Po peak. Such a component accommodates also eventual α decays in the LAr in very close vicinity to the p + detector surface. However, the number of events found with an energy >5.3 MeV is too low to fully account for the linearly modeled distribution.
The likelihood function for modeling the high-energy region dominated by α decays runs only on single-detector data, namely M1-enrBEGe and M1-enrCoax separately, in a JHEP03(2020)139 A flat prior probability is assigned to each of the fit parameters λ i . Both data sets are fit separately with a fixed bin size of 10 keV 8 as the α contamination is detector individual and the two single-detector data sets are uncorrelated in the respective energy window.
The fit results are shown in figure 9 and listed in table 8. The 210 Po component is modeled with a combination of p + contact thicknesses from 400 to 600 nm for the M1-enrBEGe data set and from 300 to 700 nm for the M1-enrCoax data set in steps of 100 nm. Further 210 Po components are rejected by a Bayes factor analysis. Impurities belonging to the 226 Ra chain are mostly located on ANG2 and thus a fit of the M1-enrCoax data set using a single p + thickness describes this component well. For the M1-enrBEGe data set we observe a very small number of counts for the 226 Ra chain, therefore, also in this case a single com- 8 The calibration curves are accurate on the sub-keV level up to the highest γ energy of about 2.6 MeV emitted by the 228 Th calibration sources. Although no major non-linearity effects were found the same accuracy cannot be guaranteed at 6 MeV. Deviations from linearity at this energy are within 10 keV, hence, we increase the bin size in the higher energy range.  ponent is sufficient. We determine a best-fit value of 100 nm and 500 nm, respectively. The estimated p-value for M1-enrBEGe is 0.2 whereas the p-value for M1-enrCoax is 0.3. The dominant spectral component below 4.5 MeV is due to degraded α events which extends down to the ROI.

B.1 Time distribution of α events
The time distribution of 210 Po decays is well known to be exponential, however, in the presence of a 210 Pb contamination a constant contribution can also be observed. 210 Pb, decaying to 210 Po, feeds a constant 210 Po component once their decay rates stabilize in a secular equilibrium. To disentangle the two we fit the time distribution of events with energies between 3.5 MeV and 5.25 MeV with a constant C and an exponential function: 2) days is the half-life of 210 Po. We use a Poisson likelihood function corrected for data acquisition dead time [50] and model the time bin content as follows  The log-likelihood can be written as a sum: We select only detectors that were ON or in anti-coincidence mode 9 in the full data taking period. In this way we avoid bias due to selection or deselection of particularly contaminated detectors. Furthermore, we exclude the initial data-taking period between December 2015 to January 2016 from the following analysis because of detector instabilities after the Phase II upgrade works. The analyzed data span from 25 th January 2016 to 3 rd April 2018 and are split into two data sets according to detector type, containing 27 enr BEGe and 7 enr Coax detectors. The fit results are shown in figure 10 and listed in

C Monte Carlo simulations and probability density functions
Background components that were identified in the energy spectra (see section 2) or in radio-purity screening measurements [5] are simulated using the MaGe software [18] based on Geant4 [19][20][21].
The Gerda Phase II detectors, their arrangement in seven strings as well as the LAr instrumentation are implemented into MaGe. A graphic rendering of the relevant implemented hardware components is presented in figure 1.
Simulations of radioactive contaminations in the following hardware components are performed: in the bulk and on the p + and n + surfaces of the germanium detectors, in the LAr, detector holder bars and plates, nylon mini-shrouds, LAr veto system (i.e. the fiber shroud, SiPMs, copper shrouds and photomultipliers) and in the signal and high-voltage flexible flat cables. The primary spectrum of the two electrons emitted in the 2νββ decay is sampled according to the distribution given in reference [22] implemented in Decay0 [23]. Note that the thickness of the detector assembly components are significantly smaller than the mean free path of the relevant simulated γ particles in the given material, thus, no significant difference can be expected between the resulting spectra of bulk and surface contaminations. The detectors n + contact thicknesses are implemented according to the values reported in references [13,14].
The 42 K decays (except for surface contaminations) are simulated homogeneously distributed in the relevant LAr volume. The following LAr volumes are chosen for the background model: the first is a cylinder centered on the detector array (h = 250 cm, r = 100 cm, simply referred to as "homogeneous" or abbreviated to "hom." in the following) subsequently divided into the volume enclosed by the mini-shrouds and the remaining one (outside the mini-shrouds); the second is a cylinder (h = 100 cm, r = 25 cm) positioned just above the array and the remaining seven are smaller cylinders (h = 20 cm, r = 5 cm), each one positioned just above each of the seven detector strings.
On top of the MaGe simulations a post-processing step is performed to compute the Probability Density Functions (PDFs) used to model the Gerda data in the statistical analysis. This includes folding in run-time dependent information, i.e. the detector status in each physics run, the finite energy resolution and threshold of each detector. All PDFs presented in the following are computed using the run-time parameters of the data sets  For the potassium tracking analysis PDFs binned in detector space are used to model the data. The rotationally symmetric single-detector PDFs for the 40 K and 42 K energy windows are shown in figure 3f and figure 12a. For two-detector events the same representation style as in figure 7 is used: projections of the two-dimensional histograms on their axis are summed, such that each two-detector event enters the final histogram twice, in the two bins associated to the respective detectors. They can be found in figure 12 together with the single-detector PDFs of the rotationally asymmetric components.
Common features can be noticed across the multitude of histogram shapes. The event rate in single-detector data is generally higher in coaxial detectors, due to their larger mass compared to BEGe detectors -maximal correlation between event rate and detector-bydetector exposure can be found in the 2νββ PDF in figure 3f. This feature is generally lost in the two-detector data: the coaxial detectors larger volume allows to stop more efficiently γ particles that would otherwise escape and eventually deposit energy in a second detector. Other similarities between different PDFs can be attributed to detectors live-times, like in the case of GD91C, which was inactive for a large fraction of the Phase II exposure and thus generally registers a low number of counts. The effects of asymmetrically distributed background contaminations are easily recognizable in the shape of the PDFs. Impurities located above the detector array are mostly seen by the upper most detectors in each string as can be seen for 40 K in the front-end electronics in figure 12a and in figure 12c and for 42 K above each mini-shroud (see figure 12e and figure 12d). Rotationally asymmetric components are mostly evident in a single string, see for example 40   (d) 40 K located close to each single mini-shroud, M2-K40 data set.

JHEP03(2020)139
All α decays in the 226 Ra to 210 Pb sub-chain and from 210 Po are simulated on the p + detector surface separately and for different thicknesses of the p + electrode. The 226 Ra chain is simulated together under the assumption that in each α decay half of the contamination is lost due to the recoil of the nucleus into the LAr. The resulting PDFs are displayed in figure 3e and figure 11a. The spectra exhibit a peak like structure with a pronounced low-energy tail. The maximum is shifted with respect to the full emission energy due to the thickness of the p + contact. The low-energy tail is characteristic for α decays; the α particle is susceptible to the change in the contact thickness when penetrating the detector surface under an incident angle and loses part of its energy before reaching the active detector volume.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.