GERDA and LEGEND: Probing the Neutrino Nature and Mass at 100 meV and beyond

The GERDA (GERmanium Detector Array) project, located at Laboratori Nazionali del Gran Sasso (LNGS), was started in 2005, a few years after the claim of evidence for the neutrinoless double beta decay (0νββ) of 76Ge to the ground state of 76Se: it is an ultra-rare process whose detection would directly establish the Majorana nature of the neutrino and provide a measurement of its mass and mass hierarchy. The aim of GERDA was to confirm or disprove the claim by an increased sensitivity experiment. After establishing the new technology of Ge detectors operated bare in liquid Argon and since 2011, GERDA efficiently collected data searching for 0νββ of 76Ge, first deploying the 76Ge-enriched detectors from two former experiments and later new detectors with enhanced signal-to-background rejection, produced from freshly 76Ge-enriched material. Since then, the GERDA setup has been upgraded twice, first in 2013–2015 and later in 2018. The period before 2013 is Phase I and that after 2015 is Phase II. Both the GERDA setup and the analysis tools evolved along the project lifetime, allowing to achieve the remarkable average energy resolution of ∼3.6 and ∼2.6 keV for Coaxial Germanium (COAX) detectors and for Broad Energy Germanium (BEGE), respectively, and the background index of 5.2+1.6 −1.3 · 10−4 cts/(keV·kg·yr) in a 230 keV net range centered at Qββ. No evidence of the 0νββ decay at Qββ = 2039.1 keV has been found, hence the limit of 1.8 · 1026 yr on the half-life (T0ν 1/2) at 90% C.L. was set with the exposure of 127.2 kg·yr. The corresponding limit range for the effective Majorana neutrino mass mee has been set to 79–180 meV. The GERDA performances in terms of background index, energy resolution and exposure are the best achieved so far by 76Ge double beta decay experiments. In Phase II, GERDA succeeded in operating in a background free regime and set a world record. In 2017, the LEGEND Collaboration was born from the merging of the GERDA and MAJORANA Collaborations and resources with the aim to further improve the GERDA sensitivity. First, the LEGEND200 project, with a mass of up to 200 kg of 76Ge-enriched detectors, aims to further improve the background index down to <0.6 · 10−3 cts/(keV·kg·yr) to explore the Inverted Hierarchy region of the neutrino mass ordering, then the LEGEND1000 (1 ton of 76Ge-enriched) will probe the Normal Hierarchy. In this paper, we describe the GERDA experiment, its evolution, the data analysis flow, a selection of its results and technological achievements, and finally the design, features and challenges of LEGEND, the GERDA prosecutor.


Introduction
In the last two decades, neutrino oscillation experiments proved irrefutably [1,2] that neutrinos oscillate into one another while propagating: this because they have a non-zero mass, and the three neutrino mass eigenstates ν i , i ∈ {1, 2, 3} with mass eigenvalues m i are a linear combination of the flavor eigenstates ν f , f ∈ {e, µ, τ}. The leading parameters driving the oscillations, i.e., the three neutrino mass differences (∆m 2 ij ) and the oscillation mixing angles θ ij , have been measured with improving accuracy [3], though not at the 1% level as in the quark sector: the Majorana CP-violating phases are still unknown, as well as the octant of θ 23 . Unfortunately, these experiments can neither establish if a neutrino is a Dirac or Majorana particle nor measure its absolute mass, and so far, have not provided information on the eigenstates mass order (hierarchy). Some (weak) preference for the normal hierarchy order, i.e., largest splitting between the second and the third mass state, has been very recently provided by SuperKamiokande [4]. In contrast, 0νββ decay experiments investigate the lepton number conservation and the nature of the neutrino (Dirac or Majorana); they measure or set limits for the absolute Majorana neutrino mass and provide information on the neutrino mass hierarchy (normal, inverted, or quasi-degenerate).
Massive Majorana neutrinos may acquire mass through the seesaw mechanism [5] and have a main role in the matter-antimatter asymmetry of the Universe.
In the Standard Model (SM), Double Beta Decay (DBD) is supposed to occur with the emission of two neutrinos (2νββ), as shown in Figure 1 (left); it was first detected on 82 Se [6] and to date has been observed in 13 nuclei [7,8]. The zero neutrino (0νββ) mode shown in Figure 1 (right), foreseen in many extensions of the SM, is so far unobserved despite the fact that it has been searched with increasing experimental sensitivity for almost 70 years: its quenching may be related to CP-violating phase cancellations, or to the quenching of the axial-vector coupling constant g A of the Nuclear Matrix Elements (NME) associated to the nuclear transition, or if neutrinos are Dirac fermions. The DBD fingerprint in experiments is a continuous spectrum ending up at the nuclear transition Q value (Q ββ ) for the 2νββ and a monochromatic line at Q ββ for the 0νββ, as shown in Figure 2: for 76 Ge the Q ββ is at 2039.15 keV. The 0νββ half-life (T 0ν 1/2 ) can be related to the effective neutrino Majorana masses (as defined by the PMNS mixing matrix U) according to [9]: (1) where M 0ν L is the NME for the light neutrino exchange (LNE) transition; M 0ν H is the NME for the heavy neutrino exchange (HNE) transition; m ββ is the effective Majorana light neutrino mass (≡ |∑ l (U 2 el m l )|); M ββ is the effective heavy neutrino mass (≡ |∑ l (U 2 el M l ) −1 |) and the factorized form of the isotopic phase space factor, G 0 ν = g 4 A G 0 ν, is used to separate the dependence on the axial-vector coupling, g A .
For each 0νββ candidate isotopes, the state-of-the-art NME calculations have a high degree of discrepancy (within a factor 2-3) largely depending on the adopted method and set of parameters, correlations functions, assumptions on intermediate states, etc. (see Table 1). High-precision experimental information from Single Charge Exchange (SCE), transfer reactions and Electron Capture (EC) as ( 3 He,t) and (d, 2 He) are used to constrain the NME calculations [10], and recently heavy ion-induced SCE and Double Charge Exchange (DCE) experiments have been proposed where the nuclear charge is changed by two units, leaving the mass number unvaried, in analogy to the ββ decay [11] nuclear model and within each model.  [20] 3.34 3. 26 2.49 In 2004, following the claim for the evidence at 4σ of 0νββ with a best value of T 0ν 1/2 = 1.19 · 10 25 yr, with a 3σ range (0.69-4.19)·10 25 yr by a part of the HDM Collaboration [21,22], the GERDA Collaboration was formed and the experimental proposal was submitted [23]: first, the design would have made it possible to reach a background index (BI) of 10 −2 cts/(keV·kg·yr) to scrutinize the controversial claim and then to reach the 10 26 yr sensitivity by further reducing the BI to 10 −3 cts/(keV·kg·yr). In the years 2005-2010, the GERDA detector and facilities were designed and built at the INFN Laboratori Nazionali del Gran Sasso (LNGS) at a depth of 3500 m w.e. (water equivalent). In those years, whereas MAJORANA further refines the background reduction techniques in the traditional approach of operating germanium detectors in vacuum, GERDA defines the technology to operate HP-Ge detectors submersed bare into liquid argon (LAr), following a suggestion by [24]; LAr serves simultaneously as a high-purity, low-Z shield against external radioactivity and as a cooling medium. Figure 3 shows the GERDA design sensitivity (S 0ν 1/2 ) on the (T 0ν 1/2 ), as a function of the exposure and for different BIs. From the exponential decay law and Poisson statistics, it directly follows that in the presence of a background, i.e., when in the experimental lifetime, a number N > 1 of events (background regime) is expected in the Region Of Interest (ROI), the T 0ν 1/2 sensitivity scales with the square root of the exposure, while when operating in a background free regime, i.e., N < 1 cts in the ROI in the experimental lifetime, the sensitivity scales linearly with the exposure: where M is the detector mass, is the estimated detection efficiency of the searched process, i.a. and A are the isotopic abundance and the atomic mass number of the exposed isotope, respectively: i.a M A gives the number of moles, and the product M t is usually defined as the exposure (E ).
The GERDA project lifetime is divided into two phases according to the background regime: in Phase I, the sensitivity was still affected by background, while Phase II was operated in a background-free regime.

The Setup
The GERDA setup [25] was designed, constructed and commissioned in the years 2005-2011 with the ambitious goal of obtaining a background index at its core of about 10 −3 cts/(keV·kg·yr). It is represented in Figure 4 where the parts introduced and/or modified in the 2015 and 2018 upgrades are labeled in red.
The core of the setup is the HP-Ge 1 detector array submersed in a 64 m 3 LAr volume complemented by a 590 m 3 purified water volume, which serves as neutron and γ absorber, contained by a stainless steel tank, on whose walls a Cerenkov muon veto detector is installed. The LAr volume of ∼75 cm diameter and ∼300 cm height surrounding the Ge array was delimited by a 30 µm thick Cu foil Oxygen Free High Conductivity (OFHC) named Rn-shroud, which prevents radons, emitted by the inner cryostat surfaces, from being transported by convection close to the Ge detector array, as shown in Figure 5.
The latter was composed of eight semi-coaxial (COAX) detectors enriched up to ∼87% in 76 Ge ( enr Ge) from the former Heidelberg-Moscow (HDM) [26] and (IGEX) [27] experiments, one coaxial nat Ge detector, and since July 2012, a string of 5 freshly produced enr Ge point-contact germanium detectors (BEGE) [28], for a total mass of ∼17.6 kg of enr Ge: two of the five BEGE assembled in a Phase I string are shown in Figure 6. The anode contacts of the COAX detectors were properly refurbished. The detector holders were made of OFHC Cu (∼80 g), and the electrical contacts (OFHC-Cu strips protected by Teflon pipe) at the detector level were spring loaded. The front-end electronics shown in Figure 6 [29] were custom developed to be low activity and operated in a cryogenic environment: each low activity PCB (see Table 2) hosting three FE circuits served three detectors. The FE PCBs were enclosed in Cu boxes ∼50 cm above the top detector of each detector string. Table 2 reports the specific radioactivity and the construction material masses in the proximity of the Ge detector array. The radioactivity has been featured by ICPMS, γ-ray spectrometry, INAA, and Rn-emanation measurements. Rn emanation of the instrumented lock and of the cryostat resulted in <10 mBq and ∼55 mBq, respectively, corresponding to a BI at the Ge-detector array of ≤ 10 −3 cts/(keV·kg·yr) on the assumption of a uniform Rn dilution in the LAr. All of the subsystems of the GERDA setup shown in Figure 4, have been fundamental for the success of the project; they served along both Phase I and Phase II. The cosmic ray panel veto, located above the clean room, integrates the water Cherenkov muon veto in the solid angle above the Ge-array in correspondence of the cryostat neck; the clean room and the lock system allowed the insertion, extraction and assembly of the Ge detector array in a safe gas Argon atmosphere, while preventing air impurities from contaminating the detectors and the LAr in the cryostat. The lock was equipped with bundles of coaxial cables woven in flat bands by Teflon wires for signal, high voltages, and power supplies; off-shell products were tested for safe operation at cryogenic temperatures and customized to minimize their radioactivity and radon emission. The cryostat and water tanks were equipped with pressure, temperature and level sensors allowing remote monitoring and safe operation of the setup. A detailed description of the GERDA Phase I setup can be found in [25].

The Early Activities and the 42 Ar Issue
Below is a selection of the early R&D and commissioning activities that paved the way for the GERDA physics achievements: • Development of the knowledge, procedures and ability to safely handle, cool down and warm up the bare HP-Ge detectors; • Findings on the role of the passivation layer for the stable operation of the bare HP-Ge detectors in LAr [31] and definition of the detector contacts (re)processing for the stable operation in LAr; • Mitigation of the 42 Ar-related 2 background: since the first commissioning, the 42 K γlines 3 intensity has been found to be higher when compared to expectations. 42 K is the decay product of 42 Ar, an isotope produced in the upper atmosphere by cosmic rays via the 40 Ar(α,2p) 42 Ar with the expected ratio N(42)/N(Ar) ∼10 −20 [32]; at the time of commissioning, the measured limit for N(42)/N(Ar) was <4.3·10 −21 [33]. GERDA found that the BI at Q ββ scaled with the 42 K γ-lines intensity. The electrostatic field dispersed in LAr from the Ge detector HV-biased surfaces was driving 42 K ions 4 close to the detector surfaces. Once in the vicinity of the detectors, the 3.5 MeV β particles travel O(1 cm) in LAr and O(2 mm) in Ge, causing background events at Q ββ mostly when entering the 1 µm thick p + contact. The 2424 and 3447 keV γ scattering in Ge generate background events at Q ββ too. To mitigate the 42 K background, the LAr in intimate contact with the Ge detectors was confined, enclosing each detector string in an OFHC Cu 60 µm thick cylinder, named "mini-shroud". It largely reduced the volume of LAr where the 42 K ions were collected and drifted to the detector anodes and to the grounded surfaces in their proximity (Cu holders). The Cu mini-shrouds allowed to reduce the BI from 18 · 10 −2 cts/(keV·kg·yr) down to 5.9 · 10 −2 cts/(keV·kg·yr) [25]. Later by improving the mini-shroud shield hermeticity and wrapping the detector HV bias contacts in Cu foils, the BI was further reduced to 2.0· 10 −2 cts/(keV·kg·yr) [25], allowing the start of physical data-taking. Concerning 39 Ar, its β endpoint at 565 keV is harmless for 0νββ decay searches, but its significant activity O(1) Bq/kg of LAr greatly reduces the GERDA potential for WIMP dark matter searches. By the spectral decomposition and fit of the physics data with the known background sources, GERDA found the 42 Ar activity to range from 72-111 µBq/kg [34] in fresh LAr distilled from the atmosphere, corresponding to a concentration N(42)/N(Ar) = 7-12· 10 −21 , exceeding the available limit [33].

Data-Taking and Treatment
The Phase I data-taking lasted from 9 November 2011 till 3 May 2013, when the exposure (E ) of 21.6 kg·yr, equivalent to 215.2 ± 7.6 moles of 76 Ge, was reached.
To compensate for the loss of one of the enr Ge-coax detectors that showed a high leakage current and to test a promising new detector type with enhanced pulse shape discrimination features [35], on July 2012, a string of enr BEGE detectors produced from a fresh batch of enr Ge material was inserted in the setup, and the data recording was restarted. Figure 7 shows the evolution of the exposure useful for the analysis, which, because of detector-related instabilities, was collected with six (of eight) enr COAX and four (of five) enr BEGE. The detector array was calibrated about every ten days, or sooner in case test pulse instabilities are observed: the black vertical lines show the physical data-taking stops for calibrations, maintenance or upgrades. For the first time in a ββ decay search experiment, a blinded analysis was performed: events falling in a 40 keV 5 region of interest (ROI) centered at Q ββ had not been reconstructed until the analysis cuts and algorithms, and their efficiencies and systematics were defined [34,36].
The data treatment is outlined in the following: first, quality cuts were applied (99.7% accepted for E> 500 keV), then single multiplicity (only one detector above the DAQ trigger) and anti-coincidence within 1 µs with the muon veto detector were required with acceptances of (94.5 ± 0.6)% and (93.7 ± 0.6)% for E > 500 keV, respectively. When considering the 100 keV energy region around Q ββ , the quality cut efficiency is unchanged within the errors, while the single-multiplicity selection keeps (66 ± 7)% of events, and the muon-veto (60 ± 7)%. This is due to the absence of the 2νββ signal, which is intrinsically single hit, and the decrease of the full containment efficiency at the increase of the γ-energy. To maximize the exposure while taking advantage of the different energy resolution, BI, and pulse shape discrimination (PSD) features of the individual detectors and time periods, data were divided into three sets: (i) GOLD COAX (E = 17.9 kg·y, BI = (1.8 ± 2) ·10 −2 cts/(keV·kg·yr)) are all the coaxial data but two runs after the insertion of BEGEs string in July 2012; (ii) SILVER COAX (E = 1.3 kg·y, BI = (6.3 +16 −14 ) · 10 −2 cts/(keV·kg·yr)) are the coaxial data from the latter two runs; (iii) BEGE (E = 2.4 kg·y, BI = (4.2 +10 −8 ) · 10 −2 cts/(keV·kg·yr)) are all the BEGE data. The evolution of the BI along the GOLDEN and SILVER data sets is shown in Figure 8. The increase of the background index in the COAX detector at the time of the pilot enr BEGE string insertion is due to a minimal quantity of Rn contamination

Filtering and Energy Calibration
The Charge Sensitive Pre-amplifier (CSP) output signals were fed to commercial digitizers based on the Analog Devices AD6645 A/D converter (14 bits, 100 MHz), and 16,384 samples were recorded. To improve the data transfer rate, for energy reconstruction only, the traces were re-binned summing up 4 consecutive bins. In this way, the waveforms contain 4096 bins of 40 ns width. After Phase I, an offline energy reconstruction was implemented; the digitized charge pulses were recorded and then analyzed with the collaboration software GELATIO [37]. The standard energy reconstruction algorithm, adopted in Phase I and for reference purposes also in Phase II, is a digital pseudo-gaussian filter. The signal was first inverted and then differentiated with a 5 µs time constant corresponding to a CR analog filter; finally, a series of 25 moving averages (MA) of 5 µs length are applied. The output signal is quasi-gaussian and its height provides the energy of the event. The 25-fold moving average, 5 µs shaping time each, maximizes the suppression of the high-frequency noise. The different steps of the filtering procedure are shown in Figure 9. Despite the fact that the pseudo-gaussian shaping is stable and relatively fast, it is limited by several factors: it can be shown that it is close to the optimal when the 1/f noise is negligible, which is not the case for GERDA, where, due to the background requirements, the CSP is located 40 to 80 cm away from the detectors. Furthermore, the adopted pseudo-gaussian shaping is the same for all detectors, although the signal profile features depend on the detector type, and the noise conditions are individual. An improvement in the energy resolution can be achieved by customizing the shaping filter parameters for each detector.
In the 0νββ experiment, a key point was to correctly address the energy scale and the resolution and monitor their stability throughout the data-taking period. The gain stability during the physical runs was monitored using a spectroscopy pulser by regular injection of charge signals to the preamplifiers input with a frequency of 0.5 Hz.
The energy calibration was performed in the case of instabilities or every 1 or 2 weeks by lowering three 228 Th sources to the vicinity of the detectors. The stability of the energy scale was monitored by comparing two consecutive calibrations: the energy shift between them was usually smaller than 1 keV at Q ββ , while the deviation of the reconstructed peak positions from the calibration curves is smaller than 0.3 keV. By the calibrations, the energy resolution is monitored as well, which was found to be also stable during the Phase I data-taking. The fit of the 1524.6 keV 42 K γ−line in the calibrated physical data is found at the proper energy but showed to be 10% broader than expected from the calibration curves. This cannot be explained by the known Doppler broadening [38] of the line and is instead related to misalignment of the calibrations of the spectra that are summed up. Taking into account the empirical 10% factor, the "detector wise-exposure weighted" energy resolution (FWHM) at Q ββ resulted in 4.8 ± 0.2 and 3.2 ± 0.2 keV for the semi-coaxial and BEGe detectors, respectively. The time profile of the FWHM at Q ββ , as measured by the 228 Th calibration series, is shown in Figure 10 for all the Phase I detectors.

Pulse Shape Analysis
The last analysis step is the event selection based on pulse shapes [36]: 1 MeV electrons travel ∼1 mm in Ge; hence, 0νββ is intrinsically a bulk, single-site event (SSE). To reject events by γ multi-site energy deposition (MSE), or by surface β and α while preserving bulk SSEs, different selection criteria and algorithms have been developed for COAX and BEGEdetectors. For COAXs, the pulse shape estimator is the output of an Artificial Neural Network (ANN), while for the BEGEs, it is the ratio of the amplitudes of the current pulse vs. charge pulse (A/E): SSE events have narrower A/E distribution than MSE ones. Events populating the 1592 keV γ-line (double escape peak of the 2614.5 keV, 208 Tl) and the full energy peaks (FEP) of the 1620 keV 212 Bi γ-line are proxies of 0νββ (SSE) and background (MSE), respectively. The acceptances are then verified on other SSE classes of events, as 2νββ and Compton edges. Figure 11 shows the distribution of the PSD estimator for different event classes for one coaxial (left) and BEGEs, respectively. When requiring an acceptance of 90% (92%) at the DEP line for the COAX (BEGEs), the acceptance of ∼50% (∼10%) for the 212 Bi γ FEP is found. The PSD systematics for the SSE selection is evaluated to be (10 ± 2)%. With the same cuts, the acceptances for the 2νββ population evaluated in the energy range from 1 to 1.45 MeV are (85 ± 2)% and (91 ± 5)% for coaxials and BEGE, respectively. The larger MSE rejection power of the BEGEs reflects the larger inhomogeneities of the carrier-driving electric field and, hence, of the drift isochrones [35]. This results in broader current pulse distribution allowing to better resolve two individual energy deposition sites.  Figure 12 shows the spectra of the three Phase I data sets, blinded at the Q ββ . 39 Ar βs dominate up to 300 keV; at this energy, in the enr Ge detectors, the 2νββ takes over up to 1.3 MeV and 40 K and 42 K γ-lines are overlapped. The prominent broad peaks from 210 Poand 226 Ra αs and their degraded event tails extending at lower energies dominate the enr COAX high energy spectrum. The enr BEGE show to be much less affected by the 210 Po contamination. γ-lines from 214 Bi and 208 Tl are present in all the three data sets. Before the unveiling, a model [34] of the radiation sources and their location was worked out. The source components, their intensity and location are identified and constrained by the characteristic γ and α-lines identified in the spectra: the model (in the following BM) reproduces the energy spectra well over almost two energy decades from 100 keV (the analysis threshold) up to 7 MeV: the minimal BM includes 2νββ in enr Ge, 39 Ar, 42 Ar and 42 K in LAr, 40 K in holders (H), responsible for both the continuum and the few visible γ lines up to energies of ∼1600 keV. 214 Bi, 228 Th, 228 Ac in detector holders, 214 Bi and 42 K βs at p + contact, 60 Co both in H and in detectors plus degraded α are the relevant components in the energy range around Q ββ and up to ∼3 MeV (Figure 13 top panel). The α region above 3 MeV is also shown in the Figure 13 bottom panel for the GOLD COAX: the data show that COAX are on average ∼10 times more contaminated in 210 Po on the p + surface and this most probably happened at the time of the COAX detector refurbishment, when the p + contact was newly implanted. The BEGE lowers contamination scales with the implanted p + surface. The BM reproduces the data when including 210 Po at the detector p+ contact, and 226 Ra both on the detector surface and in LAr. Despite the fact that BEGEs are 10 times cleaner in 210 Po than coaxials, the latter, prior PSD, have a lower background index at Q ββ (see Table 3). When more components, i.e., 228 Th from the calibration source, 214 Bi and 42 K at p + contact and 214 Bi in LAr, are included (maximal BM), the data are also well fitted, but the extra components are not strictly required.  The BM provides a solid base for the flat background hypothesis that is assumed when fitting the data.

Background Model
The unveiling confirmed that no lines are present in the data within 30 keV around Q ββ . Hence, it is correct to assume a flat background at Q ββ and estimate it by linear fitting the 230 keV energy window around Q ββ , with the exclusion of the Q ββ (2039 ± 5) keV region and of the two intervals (2104 ± 5) and (2119 ± 5) keV, corresponding to known γ-ray lines from 208 Tl and 214 Bi.

Results
The half-life on 0νββ decay is derived as with N A being Avogadro's number and m enr = 75.6 g the molar mass of the enriched material. N 0ν is the observed number of excess counts above the background or the corresponding upper limit. The efficiency , reported in Table 3, accounts for the fraction of 76 Ge atoms ( f 76 ), the active volume fraction ( f av ), the signal acceptance by PSD (ε psd ), and the efficiency for detecting the full energy peak ε f ep . The latter is the probability that a 0νββ decay taking place in the active volume of a detector releases its entire energy in it, contributing to the full energy peak at Q ββ . Energy losses are due to bremsstrahlung photons, fluorescence X-rays, or electrons escaping the detector active volume. Monte Carlo simulations yield ε f ep = 0.92 (0.90) for semi-coaxial (BEGe) detectors. Table 3. The main facts of the three Phase I data sets when (non) applying the pulse shape discrimination (PSD). "bkg" is the number of events in the 230 keV window, and BI is the respective background index, calculated as bkg/(E · 230 keV). "cts" is the observed number of events in the interval Q ββ ± 5 keV.  Table 3 reports the unveiled events (5 in the GOLDEN, 2 in the SILVER, 1 in the BEGE data sets) to be compared to the expected numbers from the flat background hypothesis, namely 5.1 in the coaxials and 2.5 in the BEGEs. The PSD rejects three events in the coaxials and a single event in the BEGEs.

Data
The results on 0νββ are obtained with PSD. A profile likelihood fit of the spectrum (shown in Figure 14) is performed; the three BIs (constants) and 1/T 0ν 1/2 (a Gaussian centered at Q ββ , σ = FWHM/2.35) are the four free parameters. The best fit returns N 0ν = 0. In a frequentist approach, this is consistent with a signal of strength N 0ν < 3.5 counts (blue line), corresponding to T 0ν 1/2 > 2.1 · 10 25 yr at 90% C.L. The claimed signal, with a halflife of T 0ν 1/2 = 1.19 · 10 25 yr [21], would produce 5.9 ± 1.4 excess counts over 2 ± 0.3 from background (red dotted line) in ±2σ around Q ββ , to be compared with the three observed. The probability of observing zero in the case of a true signal is 1%. The Bayes factor for the claimed signal (H1) versus the background only (H0) hypothesis is P(H1)/P(H0) = 2.4 · 10 −2 . Hence, GERDA does not confirm the claim, but its result alone could not rule out the full 3σ half-life range (0.6 to 4.18)·10 25 yr of Reference [21,22] and subsequent claims. When combining GERDA with HDM [26] and IGEX [27] data, the lower limit of 3.0·10 25 yr (90% C.L.) on T 0ν 1/2 was achieved and the P(H1)/P(H0) ratio of 2.0 · 10 −4 . This corresponds to m ee < 0.2-0.4 keV depending on the adopted NME 0ν and phase space factor.  [21], T 0ν 1/2 = 1.19 · 10 25 yr (red dashed) and with the 90 % upper limit derived in this work, corresponding to T 0ν 1/2 = 2.1 · 10 25 yr (blue solid). Bottom: the energy region adopted for the background interpolation.

The Upgrade of the Setup
The setup upgrade and commissioning lasted 2.5 years from summer 2013 to the end of 2015: the enr Ge mass was increased up to 35.6 kg, 20 kg being in the form of 20 BEGE detectors [28], 19 of which properly functioned, and the rest being the Phase I enr Ge-COAX detectors (15.6 kg) plus 3 nat Ge-COAX (7.6 kg). It was shown [35] that the BEGE detectors have improved pulse shape discrimination capabilities with respect to COAX. This is thanks to the highly non-uniform electric field, yielding a large spread and gradient of the holes drift velocities in the detector volume and hence longer and stretched charge integral pulses, or wider current pulses, allowing to better discriminate multi-site (MSE) energy releases. The pulse shape discrimination of Compton scattered γs was thus improved. A selection of other setup changes are listed in the following, while Table 4 reports the measured activities of the new materials and parts adopted for the Phase II setup upgrade: • The Ge-detector string was newly designed to be more compact and to minimize the assembly materials. The Ge-detector holders are now two Silicon plates connected by OFHC Cu threaded bars ( Figure 15). The detector contacts changed from spring loaded to wire bonded. Hence, each Ge detector has two evaporated Al pads to receive the wire bond; • The front-end electronics were redesigned so that one board serves four detectors [39]; this further reduces the per-channel radioactivity by a factor 1.5 and 30 for 226 Ra and 228 Th, respectively; • The whole lock cabling was renewed to reduce its Rn emanation; the adopted cables are custom made coaxials with red copper conductors and uncolored jackets, hence reducing their radioactivity and Rn emanation by a factor 15 to 25 for 226 Ra and 228 Th than in Phase I, and a factor of 50 in 40 K; • The detector contacts were made first by flex Cuflon ® and later by Pyralux ® circuits, to further reduce the mass of the Phase I contacts, and to allow wire bonding connection: these circuits showed superior performances in terms of reliability and radioactivity; • The volume of ∼50 cm diameter and ∼220 cm height was delimited by an Oxygen Free High Conductivity (OFHC) Cu foil lined with a reflector foil and equipped with 16 photomultipliers (9 top, 7 bottom): the central 100 cm of the cylinder are equipped with 800 m scintillating fibers, replacing the former Cu Rn-shroud foil, read out by Silicon Photo Multipliers (SIPM); • The Cu mini-shrouds enclosing the detector strings have been replaced by transparent mini-shroud made from the Borexino ultra-low radioactivity nylon [40] to allow the light pulse generated in the LAr in the proximity of the Ge detectors to be visible by the veto instrumentation. For this, the mini-shrouds, the fibers, the extended specular reflector (ESR) foils, and the PMTS were coated with TPB 6 . The Phase II LAr veto and PSD concepts, the array scheme, the picture of the Ge-detector array, the nylon mini-shrouds and the fiber curtain shroud, are shown in Figure 15 and described in [30].
In 2018, the setup was further upgraded (Phase II + ), replacing the nat Ge detectors with a string of enr Ge inverted coaxials (IC). Thanks to the geometry of the contacts, and hence of the charge collecting electric field, they retained the same PSD features of the BEGEs, with a significant mass of ∼2 kg and improvement of the ββ events containment efficiency and an improved signal-to-background that scales with the volume-to-surface ratio [41,42]. Moreover, to enhance the LAr-veto performances, an inner fiber shroud (shown in Figure 15) was inserted between the innermost and the surrounding Ge-detector strings. Finally, the detector Cuflon ® flex HV contacts were replaced with the Pyralux ® ones, and the signal contacts were re-routed to minimize the stray capacitance.

Data-Taking and Treatment
The Phase II data-taking time profile is shown in Figure 16; it started on 20th December 2015 and ended on the 11th November 2019, with a collected exposure of 103.7 kg·yr. The fraction of data corresponding to stable operating conditions that are used for physical analysis is about 80% of the total. The data acquisition duty cycle was ∼93% and the fraction of valid data ∼80.4%. The trigger threshold was ∼200 keV, the trigger condition being the OR of the Ge detectors. The LAr veto detector did not trigger the data acquisition, but it was read out at each Ge trigger. Gain stability and noise were monitored by test pulses injected into the front-end electronics at a rate of 0.05 Hz.
The data processing flow follows pretty much the scheme of Phase Ibut with relevant changes in each step that in cooperation with the new hardware enhanced the energy resolution, the PSD and finally, the BI. First, data are blinded in a region of ±25 keV around Q ββ . Once the energy is assigned and the energy spectra are produced (see Section 3.4), the events with E > 500 keV are processed as in Phase I [36,43]; first, quality cuts based on the flatness of the baseline, polarity and time structure of the pulse are applied to remove triggers from electrostatic micro-discharges or noise-bursts. Physical events at Q ββ are accepted with an efficiency larger than 99.9%.
Then, single detector multiplicity and anti-coincidence within 1 µs with the muon veto are required. Figure 16. The evolution of the Phase II data-taking and exposure.

Signal Denoising
To overcome the limits of the pseudo-gaussian shaping, in Phase II, an optimized shaping filter was implemented by a Finite Impulse Response (FIR) cusp-like filter with zero total area. The optimization of the energy resolution depends almost exclusively on the equivalent noise charge (ENC). Once the detector and the electronics are fixed, this is only possible by adapting the shaping filter in such a way that the ENC is minimized. It can be shown that for input referred series and parallel noise and with infinite time, the optimum shaping filter for energy estimation of a δ-like signal is an infinite cusp with the sides of the form e −t/τ s ). When dealing with waveforms of finite length and low frequency baseline fluctuations, it turns out that the best energy resolution is achieved by adopting a finite length cusp-like filter with zero total area. To overcome the cusp related ballistic deficit problem [44], which would lead to low energy tails, a flat-top has been inserted in the central part of the cusp. The resulting filter is the Zero-Area finite-length Cusp (ZAC); the construction of the filter and the final shape is illustrated in Figure 17 (left panel). The ZAC filter was then convoluted with the inverse pre-amplifier response function and this gives the filter shape shown in Figure 17 (right panel) [45]. The energy was then estimated as the maximum amplitude of the signal obtained with the convolution of the current pulse with the ZAC filter. The major improvement with respect to the pseudo-gaussian filter is that the two ZAC featuring parameters can be tailored to optimize the energy resolution against the noise condition of each detector, which may evolve over time. This is obtained offline for each calibration run and for each detector by maximizing the resolution of the 2614.5 keV γ line.

Energy Calibration and Data Partitioning
The calibrations are in fact performed by regularly exposing the HPGe detectors to three custom-made low-neutron emission 228 Th calibration sources, each of them with an activity of about 10kBq. These sources were stored within shielding above the lock system, at a vertical distance of at least 8 m to the HPGe detector array, during physical data acquisition. The sources were replaced at the beginning of Phase II to ensure a sufficient level of activity ( 228 Th half-life is about 1.9 yr). During calibration, each source was placed at three different heights to homogeneously expose the detector; the data were acquired at each location for about 30 min. Nuclei of the 228 Th isotope end up as stable 208 Pb via α and β decays with the emission of multiple mono-energetic γ rays. These result in sharp peaks in the recorded energy spectra, as shown in Figure 18, for a particular detector and calibration. The pattern of observed peaks is used to identify the known γ lines and determine both the conversion FADC to energy and resolution. From the peak positions, the calibration curve, which is the function to convert the uncalibrated energy estimator into the physical energy in keV, is obtained. To do this, the identified peak positions in terms of the uncalibrated energy estimator is plotted against their physical energies according to the literature [46] and are then fitted with a linear function. A further quadratic correction was applied to the calibration curve after the Phase II + upgrade to account for the observed deviation from a linear behavior, which is probably related to the change in cable routing. After this correction procedure, the residuals were within a few tenths of a keV. Figure 18. Combined energy spectrum for 228 Th calibration data for all enriched detectors of BEGe (in blue) and coaxial (in red) during Phase II. [47] As in Phase I, the stability of the energy scale was monitored by using both the pulser (corresponding to 3 MeV energy) and the 2614.5 keV full energy peak (FEP) in the calibration spectrum. This is to address periods with significant energy shifts. In fact, due to hardware changes or instabilities, the energy resolution and energy scale can vary over time. Therefore, the entire GERDA data period [48] was divided for each detector into stable sub-periods, hereafter referred to as partitions. The stability of each partition is determined by looking at the FWHM at the 2614.5 keV FEP to address the changes in the detector resolution and at the calibration residual at the single escape peak (SEP) to spot changes in the energy bias close to Q ββ . There are one to four partitions for each detector, although the majority of the detectors have only two partitions, the second starting at the time of the Phase II + upgrade [47]. After having combined the calibration spectra for each detector within one partition, the calibration procedure is iterated. The SEP and the DEP are excluded from the calibration line inventory, as their resolution is expected to be degraded by the Doppler effect and by events occurring in the outer regions of the detectors, respectively. Finally, the γ lines resolutions vs. the calibrated energy E is fitted by using the function σ(E) = √ a + bE. The determination of the a and b parameters allows the retrieval of the resolutions at E = Q ββ . The values shown in Figure 19 are 2.8 ± 0.3, 4.0 ± 1.3 and 2.9 ± 0.1 keV, respectively, for BEGe, coaxial and inverted coaxial detectors. The uncertainty on these values is estimated by using the standard deviation among the results obtained for the different detector partitions.  Figure 19. Effective resolution curves for BEGe (blue), coaxial (red) and IC (green) datasets [47].

Liquid Argon Veto
Next, the anti-coincidence within 6 µs with the LAr veto is required. The spectra at this level of the data treatment are shown in Figure 20. Figure 21 zooms in on the range from 500 to 1600 keV; as expected, the LAr-veto does not affect the γ-line (full energy released in the Ge detector) from the 40 K EC-decay, while it reduces ∼80% the 42 K line from the 42 K β − γ decay and the Compton edges of both lines; in fact, after the LAr-veto cut is applied, the residual spectrum < 1400 keV is perfectly fitted by the 2νββ spectrum with T 2ν 1/2 = 1.926 · 10 21 yr [49]. The triplet half-life of the LAr luminescence and the acceptance of the LAr cut along Phase II are shown in Figure 22 top and bottom panels, respectively. In 2018, at the time of the Phase II + upgrade, the cryostat volume was contaminated by an accidental small air inlet; hence, the LAr triplet half-life decreased from ∼1 µs to ∼0.9 µs. The corresponding minor increase of the physical data acceptance (from 97.7 to 98.3) may be related to this event. The acceptance values are reported in Table 5.

Pulse Shape Discrimination
Although the pulse shape analysis of GERDA Phase II follows the same approach used in Phase I (see Section 2.5), few novelties have been introduced to improve the rejection power of multi-site events. In the case of the BEGe detectors, the A/E parameter, corrected for a slight energy dependence, was used to reject background events. The low threshold for rejecting MSEs was chosen, requiring to have a 90% survival fraction of the DEP. In addition, a high threshold was defined for rejecting α events. The effectiveness of this cut in removing the αs, γ lines and Compton events is shown in Figure 23, where the physics data with and without the A/E selection are presented. The estimated 0νββ survival fraction for the BEGE Phase II data is (87.6 ± 0.1(stat) ± 2.5(syst))%. Instead, for the coaxial detectors, due to the more complicated shapes of single site events, an artificial neural network (ANN) is used for background rejection; an additional cut based on the pulse risetime (RT) is applied in order to reject the fast signals coming from surface events and due to α decays. Similarly to Phase I, the training of the PSD methods is based on the periodic 228 Th calibrations data sets. The effect of both ANN-based cut and RT cut is shown in Figure 23. Most of the α events at high energy are rejected. For the data sample relative to the coaxial detector, the estimated 0νββ survival fraction in Phase II is (84.7 ± 0.4(stat) ± 1.0(syst))%. Finally, an additional pulse shape cut was introduced for both, BEGe and coaxial detectors. The energy of background events with slow or incomplete charge collection show a bias in the energy due to the finite integration time of the ZAC filter (see Section 3.3). These events can be identified by using the difference between the reconstructed energy using the same method but different integration times, i.e., δE cut. The effect of this latest cut is shown in Figure 23 (magenta dots). The survival fraction of this cut for signal events is 99.8%. In Figure 24, the survival fractions for BEGe (upper panel) and coaxial (lower panel) detectors are shown as a function of time. The values include all the pulse shape cuts described above. Figure 23. The event distributions for the COAX and the BEGE&IC physical data sets, following the LAr veto cut (grey) and after applying the PSD classifier cut (colored) as a function of the energy: PSD for COAX is done by the ANN, the risetime and the δE selection applied one after the other. For BEGE, the PSD is performed by selecting events by the A/E and the δE, as explained in the text. The effect of the LAr cut is shown in Figure 26.

Background Model
Because of the major changes at the core of the setup, the spectral decomposition of Phase II data in featured components is newly performed [50] with a more sophisticated approach than in Phase I. The analysis is conducted prior to the application of active background suppression techniques to data, i.e., the LAr veto and pulse shape discrimination. The known background sources 42 K in LAr, 40 K 226 Ra, 238 U, 232 Th, 60 Co are located in the components within ∼80 cm from the detectors, i.e., in the detector string structural materials, the FE contacts, the mini-shrouds, LAr readout fibers and photosensors, in the FE Electronic PCB, etc. The Probability Density Functions (PDFs) used to model the different contributions are obtained from Monte Carlo simulations. The latter are performed using the MAGE simulation framework [51], based on GEANT4-v10.4 [52]. MAGE contains a software implementation of the GERDA Phase II detectors as well as the assembly and all other surrounding hardware components. Bayesian statistical analysis fits are performed to get the component intensities, and the known inventory screenings are used as priors for all background components except for α and 42 K, which are studied individually, and the results are incorporated in the full-range fit as prior distributions. The 42 K distribution inside the LAr is likely to be inhomogeneous due to drift of the ionized decay products induced by the dispersed electric field 7 . As in Phase I, the 210 Po α sources are located at the p + detector contacts. The posterior PDFs of the various components except 40 K match well with the material screening values. The runtime ON/OFF detectors and the run livetimes are taken into account, and both anti-coincidence and coincidence spectrums are simultaneously fitted. Figure 25 shows the spectral decomposition of the Phase II event energy distributions for both single-detector events (anti-coincidence spectra) and two-detector events (coincidence spectra). At low energies, the spectra are dominated by 39 Ar, 2νββ, and 40 K and 42 K γ-lines and their Comptons: at Q ββ , the background is mostly from 42 Figure 26 shows the energy distribution of GERDA Phase II events between 1.0 and 5.3 MeV, for the exposure of 103.7 kg·yr after having applied all the selection criteria. The expected distribution of 2νββ decay events with T 2ν 1/2 = (1.926 ± 0.094)10 21 yr as measured by GERDA [49] is superimposed. Figure 27 provides two zoom views around the Qββ region: in particular, the lower panel shows the 1930 to 2190 keV region where the 0νββ analysis is performed. The two shaded areas of 10 keV width around the known γ lines at 2014 keV from 208 Tl and the 2119 keV from 214 Bi are excluded from the analysis.

Results
In Figure 27, the 13 events remaining after the unblinding are shown: five are in coaxial, seven in BEGe and one in IC detectors, respectively. We point out that one event is present in the ROI (also known as "primo"), which is expected from the Poisson statistics of the measured background anyhow. The distribution of events within the analysis window is fitted, assuming a model including a Gaussian distribution for the signal, centered at Qββ having a width corresponding to the detector energy resolution plus a flat distribution for the background. The parameters retrieved from the fit procedure are the 0νββ signal strength S and the background index B. The analysis is based on an unbinned extended likelihood fit performed in both frequentist and Bayesian frameworks, as described in [48,53]. One of the improvements with respect to the analysis performed in Phase I is that the model uses different energy resolution, efficiency and exposure for each data partition (see Section 3.4), while the background and obviously the signal are common to all the partitions. The frequentist analysis gives no indication for a signal and allows setting a lower limit for the 0νββ half-life T 1/2 > 1.5 · 10 26 yr at 90% C.L. Furthermore, if the Phase I data are analyzed together with Phase II, corresponding to a total exposure of 127.2 kg·yr, the limit is T 1/2 > 1.8 · 10 26 yr at 90% C.L. The estimated value of the background from the fit is B = 5.2 +1.6 −1.3 · 10 −4 counts/(keV kg yr, which means that 0.3 counts are expected in the ROI. The analysis using the Bayesian framework has been performed by maximizing the one-dimensional posterior probability density function of the signal strength obtained by marginalizing over the other parameters. A uniform uninformative prior between 0 and 10 −24 yr −1 has been assumed for the signal S. Analyzing Phase I and Phase II data together, a limit of T 1/2 > 1.4 · 10 26 yr (90% C.I.) was found. A better limit of T 1/2 > 2.3 · 10 26 yr (90% C.I.) was obtained by using a flat prior on the Majorana neutrino mass m ee . Assuming that the decay is dominated by the exchange of light Majorana neutrinos and adopting g A = 1.27, the phase space factor and the nuclear matrix elements from literature (see Section 1), the T 0ν 1/2 limit is converted into the limit on the effective Majorana neutrino mass of m ββ < 79-180 meV at 90% C.L.

Next Generation 76 Ge Experiment: LEGEND
As a result of the successful operation and data-taking of both GERDA and MAJORANA [54], the Ge-76 neutrinoless double-beta decay community merged to build, in the next years, LEGEND [55] (Large Enriched Germanium Experiment for Neutrinoless ββ Decay). LEGEND is a phased project aiming to reach a 0νββ decay discovery potential at a half-life beyond 10 28 yr. The best technical solutions developed within GERDA and MAJORANA [56], as well as contributions from newcomers groups, will be adopted.

LEGEND-200
The first phase, i.e., LEGEND-200, will take advantage of the existing GERDA infrastructure at LNGS, which is being modified to deploy 200 kg of Ge-76 detectors in the cryostat. LEGEND-200 has a background goal of less than 0.6 cts/(FWHM·ton·yr). The achievement of this background rate will allow reaching a sensitivity above 10 27 yr with 1 ton·yr of exposure. The corresponding discovery potential is shown in Figure 28. Physical data-taking is foreseen to start by the end of 2021 to the beginning of 2022. Multiple techniques are already planned to achieve the background reduction required for LEGEND-200, such as the use of the Majorana electroformed copper, the upgrade of the GERDA LAr veto and the improvement of the front-end electronics. A crucial point will be the choice of new detectors geometry. One new option is the Inverted Coaxial Point Contact (ICPC) detector [41,42], with similar performance to the BEGe and PPC detectors and a mass as large as a coaxial detector. Five enriched ICPC detectors were produced and deployed in GERDA during May 2018. The possible design for the arrangement of the detector strings and the LAr veto in LEGEND is shown in Figure 29. A projection of the expected background contributions near Q ββ after all cuts for LEGEND-200 is shown in Figure 30 (left panel). It has been produced by using both Monte Carlo and data-driven projections of Ge U/Th, 42 K, α based on GERDA and MAJORANA data. The grey bands indicate uncertainties in assays and background rejection.

LEGEND-1000
The second stage of LEGEND will be deploying 1000 kg of Ge detectors. This increase in detector mass will require a much larger infrastructure for accommodating the cryostat vessel. An estimate of the different background contributions is necessary to reach the goal of less than 0.1 cts/(FWHM·ton·yr) for LEGEND-1000. This background reduction is required in order to achieve a 0νββ decay discovery potential at a half-life greater than 10 28 yr on a reasonable timescale (see Figure 28). The required depth to keep cosmogenic activation backgrounds (e.g., 77m Ge) within the background budget is currently under investigation and will be a contributing factor in the choice of the site where LEGEND-1000 will be located [57].

Conclusions
Fifteen years after its conception [23] GERDA, accomplished its scientific and technological goals. With an exposure of 127.2 kg·yr of HP-Ge detectors enriched up to ∼87% with the 76 Ge isotope; the world record limit on the 76 Ge T 0ν 1/2 > 1.8 · 10 26 yr has been set: when interpreting 0νββ by light Majorana neutrino (m ee ) exchange and using g A = 1.27, the limit on m ee is set at 79-180 meV depending on the adopted NME. The GERDA setup design, the novel technology of bare HP-Ge detectors operated in LAr, and the custom production of new generation enr Ge detectors with enhanced PSD features lead to top-level performances that, in synergy with the original algorithms to reject background events, have been major breakthroughs in 0νββ searches. The superior energy resolution of the Ge detectors allows for precise determination of the background sources and geometry and minimizes the leak of 2νββ events at Q ββ . Table 6 shows an up-to-date comparison of the results achieved by the world-leading projects in the field: GERDA [48] and KAMLAND-ZEN [58] set the upper limit for the T 0ν 1/2 of 76 Ge and 136 Xe at 1.8 · 10 26 yr and 1.1 · 10 26 yr, respectively; assuming the less favorable NME values, they probe the mass of the light Majorana neutrino m ee at 160-180 meV. When assuming the more favorable ones, the m ee upper limits of 60-90 meV are set by KAMLAND-ZEN [58], GERDA [48] and CUORE [59]. GERDA, with its world record background index of 5.2· 10 −3 cts/(keV·kg·yr) at the region of interest, was operated in a "background-free" regime and hence reached the nominal S 0ν 1/2 . The promising bolometric technology, with double phonon and photons read out, which is being tested by the CUPID (Cuore Upgrade with Particle Identification) projects [60,61], is pursuing the same strategy.
The knowledge of GERDA and MAJORANA on HP-Ge detector production, operation, and modelings merged into the LEGEND Collaboration, which, in 2020, took over the GERDA setup and technology and initiated the LEGEND200 and LEGEND1000 projects. At the time of this writing the LEGEND200 commissioning activities are ongoing.   Informed Consent Statement: Not applicable for studies not involving human.

Conflicts of Interest:
The authors declare no conflict of interest. 1 The Gerda HPGe detectors are made of p-type germanium. p + and n + contacts are manufactured via boron implantation and lithium diffusion, respectively. Their thickness is about 0.5-1 µm and 1 mm, respectively. The lithium-diffused n+ detector surface is a dead layer and acts as a barrier for α particles. 2 (T 1/2 = 32.9 y; β − E β = 600 keV). 3 42 K T 1/2 = 12.6 h; (β − , E β = 3.5 MeV-to 0 + of 42 Ca (81.9%), and to 2 + (17.6%), 4 + (0.05%), 3 − (0.07%), accompanied by γ emission of 1525, 2424, 3447 keV, respectively. 4 positive 42 K ions are preferentially formed as a result of the 42 Ar β decay. 5 a coarse energy calibration is provided by the FADC. 6 Tetra-Phenyl-Butadiene, an organic wavelength shifter that absorbs the 128 nm photons from the LAr scintillation light and re-emits it peaked at 420 nm with an efficiency > 95%, allows it to be collected by quartz fibers, reflected at the ESR surface and finally detected by PMTs and SiPMs.