An improved limit on the neutrinoless double-electron capture of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{36}$$\end{document}36Ar with GERDA

The GERmanium Detector Array (Gerda) experiment operated enriched high-purity germanium detectors in a liquid argon cryostat, which contains 0.33% of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{36}$$\end{document}36Ar, a candidate isotope for the two-neutrino double-electron capture (2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu $$\end{document}νECEC) and therefore for the neutrinoless double-electron capture (0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu $$\end{document}νECEC). If detected, this process would give evidence of lepton number violation and the Majorana nature of neutrinos. In the radiative 0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu $$\end{document}νECEC of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{36}$$\end{document}36Ar, a monochromatic photon is emitted with an energy of 429.88 keV, which may be detected by the Gerda germanium detectors. We searched for the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{36}$$\end{document}36Ar 0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu $$\end{document}νECEC with Gerda data, with a total live time of 4.34 year (3.08 year accumulated during Gerda Phase II and 1.26 year during Gerda Phase I). No signal was found and a 90% CL lower limit on the half-life of this process was established \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_{1/2} >1.5\cdot 10^{22} $$\end{document}T1/2>1.5·1022 year. Supplementary Information The online version contains supplementary material available at 10.1140/epjc/s10052-023-12280-6.


Introduction
The simultaneous capture of two bound atomic electrons followed by the emission of two neutrinos plus Xrays or Auger electrons, known as two-neutrino doubleelectron capture (2νECEC), is a nuclear process allowed in the Standard Model.Compared to the two-neutrino double-beta (2νββ) decay, the simultaneous emission of two electrons and two anti-neutrinos, 2νECEC processes have lower probabilities due to the smaller phase space, therefore experimentally, they are much more challenging to observe.The first direct observation of 2νECEC was made only in 2018 by the XENON1T experiment with 124 Xe [1].Previously, indications of 2νECEC were found in geochemical measurements with 130 Ba and 132 Ba [2] and in a large proportional counter experiment with 78 Kr [3].
The lepton number violating counterpart of 2νECEC, the neutrinoless double-electron capture (0νECEC), in which no neutrinos are emitted, is also predicted [4].This process must be accompanied by the emission of at least another particle to ensure energy and momentum conservation.Different modes can be considered in which 0νECEC is associated with the emission of different particles like e + e − pairs, one or two photons, or one internal conversion electron [5,6].In analogy with the neutrinoless double-beta (0νββ) decay, the 0νECEC violates the lepton number symmetry by two units and a correspondence gerda-eb@mpi-hd.mpg.deb also at: NRNU MEPhI, Moscow, Russia c now at: Duke University, Durham, NC USA d also at: Moscow Inst. of Physics and Technology, Russia e also at: Dubna State University, Dubna, Russia f also at: INFN Roma Tre, Rome, Italy implies that neutrinos have a Majorana mass component [7].Although the sensitivity of 0νECEC processes to the Majorana neutrino mass is estimated to be many orders of magnitude lower than that of the 0νββ decay, the interest in 0νECEC is theoretically motivated by the possibility of resonant enhancement when the parent nucleus and an excited state of the daughter nucleus are energetically degenerate [4,6,7,8,9].In this case, the half-life of 0νECEC processes becomes comparable to that of 0νββ decays.Experimental searches for 0νECEC have been performed by double-β decay experiments, even though with less sensitivity compared to the search for 0νββ decay [6].
The GERmanium Detector Array (Gerda) experiment, whose main goal was to search for the 0νββ decay of 76 Ge [10,11], operated enriched high purity germanium detectors in a liquid argon (LAr) cryostat, which naturally contains the 36 Ar isotope with an isotopic abundance of 0.33%. 36Ar can undergo 2νECEC to the ground state of 36 S [12].The corresponding lepton number violating process, 0νECEC, may occur via the simplest radiative mode 1   36 Ar → 36 S + γ + X K + X L . ( The 36 Ar nucleus captures one electron each from its Kand L-shells and turns into 36 S. Two X-rays are emitted, with energies E K = 2.47 keV, and E L = 0.23 keV, corresponding to the capture of the electrons from the K-and the L-shell, respectively.Given the available energy of the decay Q ECEC = (432.58±0.19) keV [14], the corresponding energy for the γ ray is E Resonance enhancement of the process is not possible for 36 Ar [6].In the light neutrino exchange scenario, assuming a Majorana mass of 0.1 eV, the half-life of 36 Ar 0νECEC is predicted in the order of 10 40 yr, with calculations based on the quasiparticle random-phase approximation (QRPA) [13].Experimental searches for 0νECEC of 36 Ar have been performed since the early stages of the Gerda experiment [15].The most stringent limit to date on the 36 Ar 0νECEC half-life is T 1/2 > 3.6 × 10 21 yr (90% C.I.), established in Phase I of the Gerda experiment [16].More recently, this process has been searched with the DEAP detector [17], although with less sensitivity than Gerda Phase I.
In this paper, we report on the search for the 429.88 keV γ line from the 36 Ar 0νECEC with the whole Gerda data, accumulated for a total live time of 3.08 yr during Gerda Phase II and 1.26 yr during Gerda Phase I. 1 Given the available energy of the process, the internal conversion mode would also be allowed for 36 Ar.Nevertheless, the latter is strongly suppressed due to argon's low atomic number and the relatively high γ energy [13].

The GERDA experiment
The Gerda experiment was located at the Laboratori Nazionali del Gran Sasso (LNGS) of INFN, in Italy [10,18,19], where a rock overburden of 3500 m water equivalent reduces the flux of cosmic muons by six orders of magnitude [10].High-purity germanium (HPGe) detectors, isotopically enriched in 76 Ge, were operated inside a 64 m 3 LAr cryostat [20].In the second phase of the experiment, 10 coaxial (including 3 detectors with natural isotopic abundance) and 30 Broad Energy Germanium (BEGe) detectors were used [18].After an upgrade in May 2018, the three natural coaxial detectors were removed, and 5 additional inverted coaxial (IC) detectors were installed [11].Detectors were mounted on 7 strings, and each string was placed inside a nylon cylinder to limit the collection of radioactive potassium ions on the detector surfaces [21].The LAr volume around the detectors was instrumented with a curtain of wavelength-shifting fibers connected to silicon photo-multipliers (SiPM) and 16 cryogenic photomultiplier tubes (PMTs) to detect scintillation light in the LAr [22,18].During the upgrade, the geometrical coverage of the fibers was improved, more SiPM channels were added, and their radiopurity increased [11].The cryostat was surrounded by a water tank containing 590 m 3 of pure water, equipped with PMTs to detect the Cherenkov light of residual cosmic muons reaching the detector site.The instrumented water tank formed, together with scintillator panels on the top of the experiment, the muon veto system [23].

Data selection
The Gerda Phase II data taking started in December 2015; it was shortly interrupted in the Summer of 2018 for the upgrade of the setup and lasted until November 2019.The total collected data used to search for the 429.88 keV γ line from the 0νECEC of 36 Ar corresponds to a live time of 3.08 yr, divided into 1.91 yr before the upgrade and 1.17 yr after the upgrade.Due to the different detector properties, e.g.energy resolution and efficiency, and the changes in the detector configuration during the upgrade, data were split into 5 data sets, namely pre-upgrade enr BEGe, preupgrade enr Coax, post-upgrade enr BEGe, post-upgrade enr Coax, and post-upgrade enr IC.The nat Coax detectors were excluded from the analysis since they have a low duty factor due to their unstable operation in Gerda Phase II and made up a minimal amount of the exposure.
Data have been processed following the procedures and digital signal processing algorithms described in [24].
The energy of an event is reconstructed using a zeroarea-cusp filter [25].Events must pass several quality cuts based on the flatness of the baseline, polarity, and time structure of the pulse to reject non-physical events.The acceptance efficiency of physical events by quality cuts is larger than 99.9% [11].Events preceded by a trigger in the muon-veto system within 10 µs are also discarded, with negligible induced dead time (<0.01%)[11].
The experimental signature used to search for 36 Ar 0νECEC in the Gerda data corresponds to the full energy deposition of the γ ray in one germanium detector.Neglecting the energy deposition of the two X-rays, no coincident energy deposition is expected, neither in the other germanium detectors nor the LAr.Consequently, the detector anti-coincidence cut and the LAr veto cut were also applied.The energy of the two X-rays is low enough that, even if they reached the germanium detector surface, they could not penetrate the 1-2 mm dead layer and, therefore, not be detected by the germanium detector.Nevertheless, since they deposit their energy in the LAr, they could be seen by the LAr instrumentation and trigger the LAr veto.The corresponding event would escape the data selection.This effect is considered in the total detection efficiency, as will be explained in section 5.The LAr veto cut reduces the background in the region of interest of this analysis by a factor of ∼ 2, as can be seen in Figure 1.In this energy region, 39 Ar β decay dominates up to the endpoint at 565 keV, while 2νββ decay is the second dominant contribution.The Pulse Shape Discrimination (PSD) cut, successfully employed in the search for 0νββ decay [26], is unsuitable for this analysis and, therefore, not used.In fact, γ rays mostly result in multiple separated energy depositions in the germanium detector, i.e. multisite events, in contrast to the single-site events produced in the 0νββ decay.In addition, the performances of the PSD cut at the energy of interest of this analysis are poorly known.Consequently, part of the data excluded in the 0νββ decay analysis from enr BEGe and enr IC data sets because of the PSD cut was instead included here.
We combine the analysis of Gerda Phase II data with that of Gerda Phase I data reported in [16].The Gerda Phase I data taking started in November 2011 and lasted until May 2013.The total collected data used for searching for 0νECEC of 36 Ar corresponded to a live time of 1.26 yr and was divided into three data sets, namely enr Coax, enr BEGe, and nat Coax.More details on the data processing and selection of these three data sets can be found in [16].It has to be noticed that the instrumentation of the LAr volume is a unique feature of Gerda Phase II and that no LAr veto cut was available in Gerda Phase I. Counts / (2 keV) 36 Ar 0νECEC 39

Energy resolution and energy scale
The energy calibration of the Gerda germanium detectors was performed during dedicated weekly calibration runs in which the germanium detectors were exposed to three 228 Th sources [27].All calibration data were combined as detailed in [27] to determine the energy scale and resolution throughout the experiment.This work uses the effective resolution curves calculated for the five analysis data sets [28].The resolution curves are evaluated at the 36 Ar 0νECEC γ energy of 429.88 keV.The energy resolution in full width at half maximum (FWHM) and their uncertainties are summarized in table 1.The uncertainty on the FWHM is calculated assuming the same relative uncertainty as for the FWHM at the Q ββ of the 76 Ge 0νββ decay (Q ββ = 2039 keV).This was calculated in [27] as exposureweighted standard deviation.The picture might be different at low energy, and the results obtained for the 0νββ decay peak at 2039 keV might not be valid for the 0νECEC peak at 429.88 keV.In fact, the lowest energy peak used to determine the resolution curves above is the 583 keV 208 Tl peak, above the energy region of interest in this analysis.To cross-check the energy resolution at the energy of interest, we use the results of the special low-energy calibration performed at the end of the Gerda data taking.This calibration run aimed to study the energy scale and stability at low energy.The energy threshold was set to 100 keV (while it was 400 keV during regular calibration runs), allowing to extend the energy range in which the resolution curve is calculated to about 238 keV, the energy of the first 212 Pb γ peak usable for the calibration.We use the peak at 583 keV as a proxy for the 0νECEC peak, being the closest in energy.We should note that also the topology of the events for the two peaks is the same.In both cases, it is a full energy deposition of the γ energy in one germanium detector, with the γ ray starting in the surrounding of the detector array.We calculate the residuals on the FWHM as the difference between the FWHM extracted in the special low-energy calibration and the value obtained evaluating the resolution curves above at 583 keV.The residuals for each detector are shown in a histogram at the left-handed side of figure 2. We find no systematic deviation of the FWHM at this energy compared to the resolution curves.The RMS of the residuals is 0.049 keV, with only one detector with a larger residual of -0.2 keV. 2he monitoring of the energy scale for the 0νββ decay search was performed using the single escape peak of 208 Tl at 2103 keV, which is typically used as a proxy for the 0νββ decay peak at Q ββ .The residuals between the peak position after energy calibration and the nominal energy value were evaluated over time, giving a mean energy bias of -0.07 keV with an average uncertainty of 0.17 keV [27].To cross-check the energy bias at the energy of interest, we use the results of the special low energy calibration run and the 583 keV peak as a proxy for the 0νECEC peak again.We calculate the residuals on the peak position as the difference between the nominal energy value and the energy value extracted from the special low-energy calibration.The residuals for each detector are shown in a histogram at the right-handed side of figure 2. We find a mean energy bias of 0.03 keV with a RMS among detectors of 0.084 keV.This is below the estimated bias uncertainty of 0.17 keV for the 0νββ decay peak at Q ββ .It should be noted that these biases are well below the binning of 1 keV used in the analysis.The effect is therefore expected to be marginal.In this work, we adopt a mean energy bias of 0 keV with an uncertainty of 0.1 keV for all the five analysis data sets.

Detection efficiency
The γ detection efficiency is defined as the probability that a 429.88 keV γ ray entirely deposits its energy inside a single germanium detector.This was determined via Monte Carlo simulations with the Geant4-based MaGe framework [29,30].In total, 10 10 γ rays with an energy of 429.88 keV were generated in a cylindrical volume of LAr, with a radius of 1.5 m and a height of 2.5 m, around the detector array.This corresponds to a net volume of LAr, after taking into account the volume occupied by the germanium detectors and structural materials, of 17.657 m 3 .The corresponding LAr mass, given the LAr density of 1385 kg/m 3 , is 24459 kg.The contribution from γ rays originating from outside this volume to the detection efficiency is negligible, as shown in figure 3. The projected distribution of vertices from which the simulated γ rays originate is shown in blue for all the events resulting in an energy deposition in the germanium detectors and black for the events resulting in the deposition of the entire 429.88 keV γ energy in one germanium detector.Only the last contribute to the γ detection efficiency, defined for each data set as the ratio between the number of events in which the full energy is deposited in one germanium detector in the specific data set and the number of initially simulated events.The number of simulated events is high enough that the statistical uncertainties on these quantities are negligible.Detector active volume and the status of each detector over the whole data taking are considered in the simulation, as detailed in [31].
The dominant systematic uncertainty on the γ detection efficiency comes from the detector active volume uncertainty.This is estimated by varying the detector dead layer in the simulation by ±1 σ, where σ is the dead layer uncertainty, and evaluating the impact on the efficiency.Typical sizes of the detector dead layers are 1-2 mm known with a typical uncertainty of 5-30 % [32].The corresponding systematic uncertainty on the γ detection efficiency is 3% for enr BEGe detectors, 4% for enr Coax detectors, and 1% for enr IC detectors.The γ detection efficiencies multiplied by the mass of LAr in the simulation volume, together with their uncertainties, are summarized in table 1 for the different data sets.The two X-rays that are emitted in the process being searched for are neglected in the simulations.As anticipated in section 3, their energy deposition in LAr could trigger the LAr veto.To account for this possibility, the survival probability of the two X-rays to the LAr veto cut is evaluated and combined with the γ detection efficiency.We use the Gerda photon detection probability map developed in [33] to estimate the probability p(x, y, z) to detect scintillation light for each simulated event starting at position (x, y, z) and corresponding to a full γ energy deposition.From this probability, the number of photons n produced by the two X-rays of total energy E X-rays = (2.47 + 0.23) keV is obtained: where 28.12 is the number of photons produced for an energy deposition of 1 keV expected in the Gerda LAr [33].The probability P that the corresponding event survives the LAr veto cut is the Poisson probability P(0, n). 3 The mean survival probability is obtained by averaging the survival probabilities of the events corresponding to a full γ energy deposition and results in P = 0.957.Thus, the data selection discards almost 5% of the events due to the X-rays depositing their energy in LAr.The calculation of the survival probability assumes that the two X-rays deposit all the energy at the exact point where the γ ray is emitted.This assumption is considered valid since the attenuation length for a 3 keV X-ray was estimated to be about 42 µm [34], negligible compared to the 3 × 3 × 3 mm 3 binning of the photon detection probability map.The main systematic uncertainty on the mean survival probability comes from the photon detection probability map.The uncertainties on this probability map given in [33] result in a 0.5% systematic uncertainty on the survival probability.Finally, we should note that the photon detection probability map assumes the pre-upgrade configuration of the LAr instrumentation [33].This means the model does not include the inner fiber shroud installed during the upgrade to improve the light detection efficiency near the germanium detectors [11].Therefore, a customized LAr veto cut was applied to select the post-upgrade data used in this work: the SiPM channels corresponding to the inner fiber shroud are not considered to build the LAr veto condition.This way, the X-rays survival probability obtained with the preupgrade photon detection probability map is extended to the post-upgrade data sets.

Analysis methods
The energy region used to set a limit on the half-life of 0νECEC of 36 Ar is defined between 410 and 450 keV (±20 keV around the γ energy of 429.88 keV, as indicated by the orange band in figure 1).Given the high statistics in this energy region, data are used in a binned form, with a 1 keV binning.It was checked that the binning choice did not impact the analysis results.In this energy region, the dominant backgrounds are the β decay of 39 Ar and the 2νββ decay of 76 Ge.Subdominant 3 Where the probability mass function for a Poisson variable is defined as contributions to the background are, in order of importance, the 42 K decays in LAr, the 40 K, 214 Pb, and 214 Bi decays in structural materials.The sum of these contributions in the analysis window can be approximated by a linear distribution, as seen in figure 1.The signal is modeled with a Gaussian peak centered at the γ energy and with the width given by the detector energy resolution (σ = FWHM/2.355).Uncertainties on the energy scale are parametrized by a shift of the signal peak δ compared to the nominal energy.
A simultaneous fit is performed on the eight data sets listed in table 1 by adopting the following binned likelihood: where the number of events in each bin is Poisson distributed, and the likelihood is given by the product of the Poisson probabilities P for all bins i and data sets d.The likelihood depends on the half-life T 1/2 of the investigated process, which is a common parameter among the eight data sets and is the only parameter of interest, and on some nuisance parameters ⃗ θ that are data set specific and affect both the signal and background distributions.Gaussian pull terms Pull( ⃗ θ) are introduced in the likelihood to constrain some of the nuisance parameters.Finally, n di denotes the number of observed events in the data set d and bin i, and µ di is the expectation value for the same data set and bin.The latter is given by the sum of the signal and background in that bin: µ di = b di + s di .The number of signal events s di is given by the integral of the signal distribution for the data set d in the bin i.This is a Gaussian distribution centered at E + δ d (E), where E is the γ energy of 429.88 keV and δ d (E) the energy bias for the data set d calculated for the same energy, and with the width given by the detector energy resolution σ d (E) = FWHM(E)/2.355evaluated for the same data set and at the same energy.The total number of signal events in a data set d is related to the half-life T 1/2 through the relation: where N A is the Avogadro constant, M 36 is the molar mass of argon (35.968 g/mol), m LAr,d is the mass of LAr in the simulations from which the γ detection efficiencies are extracted (the product ε γ • m LAr is given in table 1 for each analysis data set), f 36 is the abundance of 36 Ar in ultra-pure natural Argon (0.334%) [35], and t d is the live time of the experiment.The total efficiency ε tot,d for the Phase II data sets is given by the product , where ε γ is the γ detection efficiency, ε X the X-rays survival probability (both discussed in section 5), and ε LAr is the efficiency of the LAr veto cut.The latter was estimated to be (97.7±0.1)% for the pre-upgrade data and (98.2±0.1)% for the post-upgrade data [11].The total efficiency of the Phase I data sets equals the γ detection efficiency ε γ , because no LAr veto cut was available in Gerda Phase I. Analogously, the number of background events b di is given by the integral of the background distribution for the data set d in the bin i.The background distribution is a linear function that depends on two parameters, the normalization and the slope, both data set-specific.We verified that the first-order polynomial function describes the data in this energy region well and that a second-order polynomial function does not fit the data better.In modeling the background of Phase I data, an additional Gaussian distribution is used to describe the full energy deposition of the 433.9 keV γ ray from 108m Ag, which lies in the energy region of the analysis.Contamination from 108m Ag was observed in the screening measurements, and all the three expected γ lines from 108m Ag were observed in Gerda Phase I data [16,36].The origin of the 108m Ag contamination in Gerda Phase I was found in the signal cables [36], which were exchanged in Gerda Phase II [18].In addition, none of these γ lines was observed in Gerda Phase II data after the LAr veto cut.The decay of 108m Ag proceeds through a cascade of three equally probably γ rays at energies of 433.9 keV, 614.3 keV, and 722.9 keV.Therefore, even if any 108m Ag contamination were still present in Gerda Phase II, the LAr veto cut would likely discard the corresponding events.In total, the fit has 42 floating parameters, 22 describing the signal peak (ε d , δ d , σ d ) 4 , 10 for the linear background of Phase II data sets, 6 for the linear background of Phase I data sets, 3 parameters for the number of 108m Ag events in Phase I data sets, plus one common parameter to all data sets T −1 1/2 .The latter is constrained to positive values.Gaussian pull terms in the likelihood given in Eq. 3 constrain some of the nuisance parameters, namely the efficiency ε d , the energy bias δ d , and the energy resolution σ d around their central value and uncertainty.All the other nuisance parameters are free and unconstrained, and their uncertainties are propagated into the result by profiling.
To set a lower limit on the half-life of the investigated process, we use a modified frequentist approach, namely the CL s method [37].The latter was found to be a more appropriate choice in the case of an experiment with low sensitivity or, in different words, a background-dominated experiment [37].Compared to a pure frequentist approach, the CL s exclusion region does not assure the correct coverage and often results in an over-coverage, thus a more conservative result.The profile likelihood ratio test statistic is used for the p-value calculation.Asymptotic distributions of the test statistic and the Asimov data set are used [38].The statistics in each bin is high enough for this assumption to be valid.

Results
The best fit, defined as the minimum of the profiled likelihood ratio, yields T −1 1/2 = 0, i.e. we do not observe any signal events from 0νECEC.Data from the five Gerda Phase II analysis data sets and in the energy region of the analysis are shown in figure 4  Systematic uncertainties on the efficiency ε d , the energy bias δ d , and the energy resolution σ d are identified as primary sources of systematic uncertainties and included in the likelihood through nuisance parameters constrained by Gaussian pull terms as explained in section 6.Their overall effect on the limit derived in section 7 is estimated to be 2%.Potential systematic uncertainties related to the fit model, particularly the background distribution, are also investigated.First, the assumption of a linear distribution is compared to a more general second-order polynomial distribution.This has a negligible impact on the result.The presence of additional structures in the background is also investigated.As discussed in section 6, a γ line from 108m Ag, very close to the expected signal energy, is included in the background model of the Phase I data sets, as in previous analysis [16].A possible systematic uncertainty due to the above γ line in Phase II data is investigated by introducing it in the background model.This would worsen our result of a 2%.

Conclusions
In this work, we searched for the 429.88 keV γ line from the 36 Ar 0νECEC using the final total exposure of the Gerda Phase II experiment, combined with the Gerda Phase I exposure.No signal was observed, and a lower limit on the half-life of this process was derived, yielding T 1/2 > 1.5 • 10 22 yr (90% C.L.).This is the most stringent limit on the half-life of the 36 Ar 0νECEC.This work shows that the potential of the Gerda experiment in investigating physics beyond the Standard Model extends further than the search for the 0νββ decay of 76 Ge (see also [39,40]).Even if the sensitivity is many orders of magnitude below the theoretical expectation for this process, to our knowledge, the Gerda experiment was, to date, the only experiment with the capability to search for the 0νECEC of 36 Ar with competitive sensitivities.The Gerda sensitivity is limited by the physical background from 39 Ar β and 76 Ge 2νββ decays in the energy region where the γ peak is expected, which is, for instance, orders of magnitude higher than the background in the region of interest for the 76 Ge 0νββ decay.An additional limiting factor is the low detection efficiency since the γ ray is emitted in the LAr and must be detected in one of the germanium detectors.Only γ rays emitted in the proximity of the detector array contribute to the total efficiency as discussed in section 5 (See figure 3).
Among the planned future experiments, the Large Enriched Germanium Experiment for Neutrinoless-ββ Decay (LEGEND) experiment can extend the search for the 0νECEC of 36 Ar to higher sensitivity.In the first phase of the project, LEGEND-200 will deploy about 200 kg of germanium detectors.This is more than a factor of four compared to the Gerda detector mass and will imply a higher detection efficiency to the γ ray emitted in this process.On the other hand, the background in the energy region where the γ peak is expected should be comparable to the Gerda back- Observed CLs Expected CLs 68% interval 95% interval Fig. 5 CLs as a function of the inverse of the half-life obtained in the analysis of the combined Gerda Phase I and Phase II data.The median of the CLs distribution for the Gerda experiment under the no signal hypothesis and the observed CLs for the Gerda data are shown by the continuous black line and the dashed line, respectively.The spread of the CLs expected distribution, given by the 68% and 95% probability intervals, is also shown by the colored bands.The 90% C.L. limit (sensitivity) is given by the solid (dashed) black line intersection with the dotted line, corresponding to a CLs of 0.1.
ground, largely dominated by the 39 Ar β decay.Still, an improvement in the current sensitivity is foreseen.LEGEND-1000 will deploy about 1 ton of germanium detectors, implying an even higher detection efficiency to the γ ray emitted in this process.In addition, using underground Ar instead of atmospheric Ar is intended.This is depleted of 39 Ar, which is the main background contribution in this search.A significant improvement in the sensitivity is therefore expected.To our knowledge, no other planned experiment has competitive sensitivity to LEGEND in the search for 0νECEC of 36 Ar.

Fig. 2 (
Fig. 2 (Left) Distribution of the energy resolution (FWHM) residuals for the 583 keV calibration peak.(Right) Distribution of the peak position residuals for the same calibration peak.The mean and the RMS of the two distributions are indicated.

Fig. 3
Fig.3Projected distribution of vertices from which the simulated γ rays originate.γ rays with an energy of 429.88 keV are simulated uniformly in the cylindrical volume.Only those originating from the blue vertices deposit some energy in the germanium detectors, while only those originating from the back vertices deposit the entire energy in one germanium detector, thus contributing to the γ detection efficiency.
together with the best-fit model and the residuals normalized to the expected statistical fluctuations of the bins.The 90% C.L. limit on the half-life is obtained by scanning the observed CL s over different values of T −1 1/2 and finding the value for which CL s = 0.1.For Gerda Phase II data only, this gives T 1/2 > 1.3 • 10 22 yr.The 90% C.L. sensitivity of the Gerda Phase II experiment, i.e. the median expectation under the no signal hypothesis, is obtained analogously by scanning the expected CL s over different values of T −1 1/2 and finding the value for which CL s = 0.1.The latter gives T 1/2 > 8.0 • 10 21 yr.The analysis of the combined Gerda Phase I and Phase II data gives a 90% C.L. sensitivity of T 1/2 > 8.6 • 10 21 yr and an observed lower limit of T 1/2 > 1.5 • 10 22 yr. Figure 5 shows the scan of the observed and expected CL s over a range of values of T −1 1/2 obtained in the analysis of the combined Gerda Phase I and Phase II data.

3 Fig. 4
Fig.4 Best fit of the combined Gerda data.The blue line shows the combined best-fit model, corresponding to T −1 1/2 = 0.The dashed orange line indicates the energy at which a γ line from 0νECEC is expected, and the orange peak displays the expected signal for a half-life equal to the 90% C.L. lower limit 1.5 • 10 22 yr.The pulls, i.e. residuals normalized to the expected statistical fluctuations of the bins, are shown in the bottom panels for each data set.
Energy distribution of the low energy Gerda Phase II data before and after LAr veto cut.The left part of the spectrum is dominated by the 39 Ar β decay with an endpoint at 565 keV.On the right side, the 2νββ decay dominates.Some known γ lines are visible and labeled.The orange dotted line indicates the energy at which the 36 Ar 0νECEC is expected, and the orange band indicates the energy region used in the analysis.

Table 1
[16]gy resolution (FWHM) and γ detection efficiency (multiplied by the simulated mass of LAr) for the analysis data sets.The values for the Phase I data sets are taken from[16].