Results on $\beta\beta$ decay with emission of two neutrinos or Majorons in $^{76}$Ge from GERDA Phase I

A search for neutrinoless $\beta\beta$ decay processes accompanied with Majoron emission has been performed using data collected during Phase I of the GERmanium Detector Array (GERDA) experiment at the Laboratori Nazionali del Gran Sasso of INFN (Italy). Processes with spectral indices n = 1, 2, 3, 7 were searched for. No signals were found and lower limits of the order of 10$^{23}$ yr on their half-lives were derived, yielding substantially improved results compared to previous experiments with $^{76}$Ge. A new result for the half-life of the neutrino-accompanied $\beta\beta$ decay of $^{76}$Ge with significantly reduced uncertainties is also given, resulting in $T^{2\nu}_{1/2} = (1.926 \pm 0.095)\cdot10^{21}$ yr.


Introduction
Neutrinoless double beta (0νββ) decay is regarded as the gold-plated process for probing the fundamental character of neutrinos. Observation of this process would imply total lepton number violation by two units and that neutrinos have a Majorana mass component. Although the main focus of experimental efforts lies on the detection of 0νββ decay mediated by light Majorana neutrino exchange, there are also many other proposed mechanisms which are being searched for. Some exotic models predict 0νββ decays proceeding through the emission of a massless Goldstone boson, called Majoron. Predictions of different models depend on its transformation properties under weak isospin, singlet [1], doublet [2] and triplet [3]. Precise measurements of the invisible width of the Z boson at LEP [4] greatly disfavour triplet and pure doublet models. Several new Majoron models have been developed subsequently in which the Majoron carries leptonic charge and cannot be a Goldstone boson [5,6] or in which the 0νββ decay proceeds through the emission of two Majorons [7].
All these models predict different shapes of the two emitted electrons' summed energy spectrum. The predicted spectral shapes are essentially defined by the phase space of the emitted particles: a also at: Moscow Inst. of Physics and Technology, Russia b also at: Int. Univ. for Nature, Society and Man "Dubna", Dubna, Russia c Correspondence, email: gerda-eb@mpi-hd.mpg.de where K is the summed energy of the two electrons, G is the phase space, Q ββ is the Q value of the 0νββ decay and n is the spectral index of the model. Single Majoron emitting ββ decays can be roughly divided into three classes, n = 1, n = 2, and n = 3. Double Majoron emitting decays can have either n = 3 or n = 7. Their characteristic spectral shapes differ from that of twoneutrino ββ decay (2νββ), for which n = 5. This allows for discrimination between the processes. Experimental searches for ββ decay mediated by emission of one or two Majorons (0νββχ) have been performed by the Heidelberg-Moscow experiment (HdM) for 76 Ge [8,9]; by Nemo-2 and Nemo-3 for 100 Mo, 116 Cd, 82 Se, 96 Zr, 130 Te [10,11]; by ELEGANT V for 100 Mo [12]; by DAMA [13] and by KAMland-Zen for 136 Xe [14]. None of these experiments have seen an excess of events that could be interpreted as a Majoron signal; they reported lower limits on the half-lives of the processes that involve Majoron emission.
The 2νββ decay process conserves lepton number and is independent of the nature of the neutrino. It has been detected for eleven nuclides so far, with measured half-lives (T 2ν 1/2 ) in the range of 7×10 18 −2×10 24 yr [15,16,17]. The knowledge of T 2ν 1/2 allows for extraction of the nuclear matrix element, M 2ν , which can provide some constraints on that of 0νββ decay, M 0ν , if the evaluations of M for the two processes are performed within the same model [18,19].
This paper reports on the search for neutrinoless double beta decay of 76 Ge with Majoron emission (0νββχ) and a new analysis of the half-life of the 2νββ decay of 76 Ge using data collected by the Gerda experiment during its Phase I. 2νββ decay is a well established and previously observed process, while 0νββχ decay is a hypothetical one. In the first case the half-life is extracted, while for the second one a limit is set. This leads to slightly different approaches in the analyses leading to different data sets and background components being used.

The GERDA experiment
The main aim of the Gerda experiment [20] at the Laboratori Nazionali del Gran Sasso (LNGS) of INFN in Italy is to search for 0νββ decay of 76 Ge. The core of the setup is an array of high-purity germanium (HPGe) detectors made from isotopically modified material with 76 Ge enriched to ∼86 % ( enr Ge), mounted in low-mass copper supports (holders) and immersed in a 64 m 3 cryostat filled with liquid argon (LAr). The LAr serves as cooling medium and shield against external backgrounds. The shielding is complemented by water in a tank of 10 m in diameter which is instrumented with photomultipliers to detect Cherenkov light generated in muon-induced showers [20].
The array of HPGe detectors is arranged in strings. Each string is enclosed with a cylinder, made from 60 µm thick Cu foil, called mini-shroud, to mitigate the background coming from the decay of 42 Ar present in the LAr. Moreover, in order to prevent contamination from radon within the cryostat, a cylinder, made from 30 µm thick Cu foil, called radon-shroud, separates the central part of the cryostat, where the detectors are located, from the rest. The HPGe detector signals are read out with custom-made charge sensitive preamplifiers optimized for low radioactivity, which are operated close to the detectors in the LAr. The analog signals are digitized with 100 MHz Flash ADCs (FADC) and analyzed offline. If one of the detectors has an energy deposition above the trigger threshold (40-100 keV), all channels are read out. Reprocessed p-type coaxial detectors from the HdM [21] and Igex [22] experiments were operated together with Broad Energy Germanium (BEGe) type detectors manufactured by Canberra [23,24].
As explained in section 5, some background components have different effects on the two detector types due to their peculiar geometry. A schematic drawing of a coaxial detector type is shown in the top part of Fig. 1 [25,26]. For the coaxial detectors, a data set collected for 1.3 kg·yr exposure during a restricted time period around the deployment of the BEGe detectors is discarded due to a higher background level. Also one of the coaxial detectors, RG2, is not considered for the data analysis starting from March 2013, as its high voltage had to be reduced below depletion voltage due to increased leakage current. The energy calibration of the detectors was performed using the information from dedicated calibration runs. For these calibration runs, three 228 Th sources were lowered to the vicinity of the detectors. The stability of the energy scale was monitored by performing such calibration runs every one or two weeks. Moreover, the stability of the system was continuously monitored by injecting charge pulses into the test input of the preamplifiers. Using physics data, the interpolated FWHM values at Q ββ averaged with the exposure are (4.8 ± 0.2) keV for the coaxial detectors and (3.2 ± 0.2) keV for the BEGe detectors.
All steps of the offline processing of the Gerda data were performed within the software framework Gelatio [27]. The energy deposited in each detector was extracted from the respective charge pulse by applying a approximate Gaussian filter [28]. Non-physical events, such as discharges, cross-talk and pick-up noise events, were rejected by quality cuts based on the time position of the rising edge, the information from the Gaussian filter, the rise time and the charge pulse height, which must not exceed the dynamic range of the FADCs. Pileup and accidental coincidences were removed from the data set using cuts based on the baseline slope, the number of triggers and the position of the rising edge. The rate of pile-up and accidental coincidence events is negligible in the Gerda data due to the extremely low event rate. The loss due to mis-classification by the quality cuts was <0.1 % for events with energies above 1 MeV. All events that come within 8 µs of a signal from the muon veto were rejected. Finally, only events that survive the detector anti-coincidence cut were considered. This means, that all events with an energy deposition > 50 keV in more than one detector in the array were not taken into account. Since 2νββ and 0νββχ events release their energy within a small volume inside the detectors, almost no signal events were lost by this cut, while a part of the γ-induced background events were rejected.

Analysis Strategy
The two analyses described in this paper are different in the sense that for 2νββ decay a parameter is extracted for a well established and known process, while in the case of the search for 0νββχ decay limits for a hypothetical process are set. In order to minimize the systematic uncertainties for the extraction of the T 2ν 1/2 it is favorable to use a well defined and controlled subset of the data and to use only well identified background processes. For 0νββχ limit setting it is favorable to maximize the exposure and to take into account all known possible background processes that can not be unambiguously detected but could mimic 0νββχ decay.
For the T 2ν 1/2 analysis the golden data set (17.9 kg·yr) with the coaxial detectors is used in order to have a large data sample obtained in well controlled experimental conditions. The Majoron analysis uses both the golden data set and the BEGe data set for a total exposure of 20.3 kg·yr in order to maximize the sensitivity.
The background model for the T 2ν 1/2 analysis uses a minimal number of components, assuming all sources near to the detectors [25,29]. For the Majoron analysis, an expanded model is used [30], taking into account also additional medium and far distant positions for some of the sources. This becomes necessary when searching for rare processes such as Majoron emission, where all possible sources of background which could simulate the exotic process have to be considered. Therefore, even the slight differences resulting, for example from a variation of the source position, have to be evaluated.
In both analyses, the experimental spectra of the coaxial and BEGe detectors are analyzed using the Bayesian Analysis Toolkit (Bat) [31].

The background model
The background sources considered in the models were identified by their prominent structures in the energy spectra and were also expected on the basis of material screening measurements. The spectral shapes of individual background contributions were obtained by using a detailed implementation of the experimental setup in the Monte Carlo (MC) simulation framework MaGe [32]. A Bayesian spectral fit of the measured energy spectrum with the simulated spectra was performed in an energy range from 570 keV up to the end of the dynamic range at 7500 keV. The low energy limit is motivated by the β-decay of 39 Ar, which gives a large contribution up to its Q β -value of 565 keV.
The following background components were used for the extraction of the T 2ν 1/2 (minimum model in Refs. [25,29]): (1) 76 Ge 2νββ decay, (2) 214 Bi, 228 Ac, 228 Th, 60 Co and 40 K decays in the close vicinity of the detectors (<2 cm, represented by decays in the detector holders in the MC simulation), (3) decays of 60 Co inside the detectors, constrained by the maximum expected activity from their cosmogenic activation history, (4) 42 K decays in LAr assuming a uniform distribution, (5) αmodel that accounts for α decays originating from 210 Po and 226 Ra contaminations on the p + surface of the detectors as well as from 222 Rn in the LAr, and finally (6) 214 Bi decays on the p + surface, constrained by the estimated 226 Ra activity from the α-model.
The parameters of all components besides the constrained ones were given a flat prior probability distribution. There are no strong correlations between the model parameters since all considered background components have characteristic features such as γ-ray lines or peak-like structures at different energies. The ratios of the γ-ray line intensities from the individual considered background sources suggest contaminations dominantly in locations close to the detectors. Hence, the minimum model takes into account only the close-by source locations. Nevertheless, the screening measurements indicate contaminations of materials in farther locations as well. An additional contribution can come from 42 K decays at or near the detector n + surfaces (see Fig. 1) with a specific activity higher than that for the uniform distribution assumption. This component is the dominating one for the BEGe data set, as the thinner dead layer thickness of BEGes of roughly 1 mm allows penetration of the electrons emitted in the decay of 42 K to the active volume, while for coaxial detectors the dead layer thickness of ∼2 mm efficiently shields this background component.
The spectral shapes of the contributions from the background sources without significant multiple γ peaks at different source locations differ only marginally. This makes it impossible to pinpoint the exact source locations given the available statistics of the measured spectra. Therefore, variations of the source locations for the considered decays were taken into account when evaluating the systematic uncertainty on T 2ν 1/2 . For the Majoron analysis additional background components were used [30], including also medium and far distant contributions. For the coaxial detectors 42 K on the n + and on the p + contacts was added to the list of the close sources of the previous background model. For medium distances, i.e. between 2 cm and 50 cm from the detectors, contributions from the following sources were added: 214 Bi, 228 Th and 228 Ac. A 228 Th contamination was chosen as a representative for far distant sources (above 50 cm). Whenever possible, screening measurements were used to constrain the lower limit of the expected background events.
In the Majoron analysis, also the data collected with the BEGe diodes were used in order to maximize the exposure. Consequently, the background model developed for these detectors was used [25,30]. The same close, medium and far distant sources as for the coaxial detectors were used. 68 Ge was added as internal source. This was necessary in order to take into account the cosmic activation of the germanium due to the recent production of these diodes.
6 Determination of the half-life of 2νββ decay

Analysis
The T 2ν 1/2 of 2νββ decay of 76 Ge was determined considering the golden data set of Phase I, amounting to an exposure of 17.9 kg·yr, and using the background model prediction for the contribution of the 2νββ spectrum to the overall energy spectrum. Details of the background analysis can be found in Ref. [29].
The global fit for the background modeling was performed on the summed energy spectrum of the coaxial detectors using a bin width of 30 keV. Thus, the scaling parameter of the 2νββ spectrum in the model, N fit 2ν , gives the number of events in the 2νββ spectrum in the fit window of 570-7500 keV for all detectors. Using this result for the number of measured 2νββ events, the half-life is calculated as where N A is Avogadro's constant and m enr = 75.6 g is the molar mass of the enriched material. The summation runs over all the detectors (N det ) considered in the data set. All detector related parameters like the detector mass (M i ), the time of the data taking for each detector (t i ), the fraction of 76 Ge atoms (f 76,i ), the active volume fraction (f AV,i ), and the detection efficiencies in the active volume (ε fit AV,i ) and in the dead layer (ε fit DL,i ) are taken into account separately for the individual detectors. All values are listed in Table 1. The efficiency ε fit AV,i (ε fit DL,i ) corresponds to the probability that a 2νββ decay taking place in the active volume (dead layer) of the detector deposits detectable energy in the fit window considered for the background model. The detection efficiencies, on average ε fit AV = 0.667 and ε fit DL = 0.011, are obtained through dedicated MC simulations. The statistical uncertainty due to the number of simulated events is on the order of 0.1 %.
The background model resulted in a scaling parameter of N fit 2ν = 25690 +310 −330 for the 2νββ spectrum, which is the best fit parameter. The uncertainty is given by the smallest 68 % probability interval of the marginalized posterior probability distribution. Using this result, the half-life derived according to Eq. 2 is

Systematic Uncertainties
The systematic uncertainties affecting the results for T 2ν 1/2 were grouped into the three categories (i) detector parameters and fit model, (ii) MC simulation, and (iii) data acquisition and selection. The contributions to the total systematic uncertainty on T 2ν 1/2 are summarized in Table 2.
(i) detector parameters and fit model -The systematic uncertainty on the active 76 Ge exposure (E AV,76 ) was determined using a MC approach. E AV,76 is defined as For evaluating its uncertainty, the parameters of the individual detectors were randomly sampled from Gaussian distributions with mean values and standard deviations according to the corresponding values listed in Table 1. The correlated terms for f AV were also taken into account. The uncertainty on the live time t is 0.3 %, whereas the total detector masses are known with good accuracy (uncertainty smaller than 0.1 %). The calculation yields E AV,76 = (13.45 ± 0.54) kg·yr. The uncertainty of 4 % is driven by the uncertainties on f AV and f 76 , which mainly affect the number of 76 Ge nuclei in the active volume of the detectors, with a relatively smaller impact on the detection efficiency for the background sources. -The reference background model used for determining T 2ν 1/2 accounts only for the dominant source locations in the setup. The systematic uncertainty due to the choice of the background model components was evaluated by repeating the global fit with alternative models, which account for different source locations for all the background sources considered in the reference model. The model that accounts for 228 Th and 228 Ac contributions also in the radonshroud instead of only in the holders results in a 1.4 % longer T 2ν 1/2 . The same increase occurs if 40 K in the radon-shroud is added to the model components. The model including the contribution from 214 Bi in the radon-shroud in addition to the p + surface and holders yields a 0.7 % longer T 2ν 1/2 . In all the cases mentioned above, the contribution from background in the 2νββ spectrum region increases, since the peak-to-Compton ratio of the γ-rays decreases for farther source locations leading to longer T 2ν 1/2 estimates. Excluding contributions from very close source locations, like 214 Bi on the p + surface and 60 Co on the germanium, results in a smaller increase of the best T 2ν 1/2 estimate. In this case, the contributions from these components are compensated by 214 Bi and 60 Co decays in the holders, respectively. Consequently, the source locations are moved further out with respect to the reference model. Consistently, the models that include additional contributions from close source locations yield a decrease in the T 2ν 1/2 value, e.g. including 214 Bi in LAr close to the p + surface (-1.0 %) or 42 K on the n + (-1.2 %) and p + (-0.6 %) surfaces. Comparing alternative background models to the reference one, the deviations in the T 2ν 1/2 result range between -1.2 % and +1.4 %.
-For the standard fit, a bin width of 30 keV was used for the data and MC energy spectra. In order to take into account the systematic uncertainty related to binning effects, the fit was repeated twice using bin widths of 10 and 50 keV. The bin width of 10 keV was chosen in order to minimize as much as possible the bin size taking into account the energy resolution of ≈4.5 keV of the coaxial detectors and the necessity to have enough statistics in all bins. Above 50 keV, peak structures are washed out, leading to a deterioration of the fit. The deviations in the T 2ν 1/2 result range between -0.5 % and +0.5 % with respect to that using the standard bin width. -The primary spectrum of the two electrons emitted in the 2νββ decay of 76 Ge, which was then fed into the MC simulation, was sampled according to the distribution given in Ref. [33] implemented in De-cay0 [34]. The systematic uncertainty due to the assumed 2νββ spectral shape was evaluated by comparing the spectrum generated by Decay0 to the one given in Ref. [35]. Considering the analysis window used for background modeling, the maximum deviation is 0.2 % and the total deviation of the integral in the analysis window is 0.1 %. When the fit with the background model is repeated using the spectrum of Ref. [35], the difference from the reference T 2ν 1/2 result is less than 0.1 %. -A possible effect of a transition layer, where it is assumed that the n + dead layer on the detector surfaces is partially active, has been investigated [36,37]. The dead layer thickness for individual detectors assumed in MC simulations were given according to the listed values in Ref. [25]. The transition layer is modeled using two different assumptions: a linearly and an exponentially increasing charge collection efficiency in the dead layer. The systematic uncertainty on T 2ν 1/2 due to the 2νββ spectrum simulated with the transition layer is found to be negligible.

(ii) MC simulation
The uncertainty related to the MC simulation arises from the precision of the experimental geometry model implemented in MaGe (1 %) and from the accuracy of particle tracking (2 %) performed by Geant4 [38,39]. The total MC simulation uncertainty was estimated to be 2.2 % by summing in quadrature the aforementioned contributions.

(iii) Data acquisition and selection
The trigger and reconstruction efficiencies for physical events are practically 100 % above 100 keV in Gerda. The performance of the quality cuts applied in Phase I data has been investigated through a visual analysis. The total uncertainty related to data acquisition and selection was estimated to be less than 0.1 %.
Summing in quadrature the uncertainties of the three groups gives a total systematic uncertainty of ±4.8 %.

Results and Discussion
with the latter combining in quadrature the statistical (fit) and systematic uncertainties. The total uncertainty of 4.9 % is dominated by the systematic uncertainties.
The largest contribution to the systematic uncertainties comes from the uncertainty on the active 76 Ge exposure (4 %), which can only be reduced by performing new and more precise measurements of the active masses of the coaxial detectors. Other significant contributions are related to the Monte Carlo simulations (2.2 %) and to the background model assumptions ( +1.4 % −1.2 % ). The latter have been significantly reduced in this analysis compared to the analysis of the first 5 kg·yr of Phase I data reported in Ref. [41], where the systematic uncertainty due to the background model was +5.7 % −2.1 % . The new result is in good agreement with that mentioned above. Adding further identified components to the reference background model results in a slight increase of the best T 2ν 1/2 estimate. The background level achieved in Gerda Phase I is about one order of magnitude lower with respect to predecessor 76 Ge experiments, and has allowed the measurement of T 2ν 1/2 with an unprecedented signal-tobackground ratio of 3:1 in the 570-2039 keV interval. The ratio amounts to 4:1 for the smaller interval of 600-1800 keV.

Analysis
The search for 0νββχ was performed using the golden and BEGe data sets, amounting to a total exposure of 20.3 kg·yr. The analysis employed the background model described in section 5. The information from the two data sets was combined in one fit, while keeping their energy spectra distinct. A separate fit was performed for each spectral index, containing the background contributions, the contributions from 2νββ decay, and also the Majoron component under study. A single parameter, T 0νχ 1/2 , is considered common for the two data sets. It is defined as the half-life of the respective Majoron accompanied mode.
In order to improve the detection efficiency for the Majoron processes with low n (n = 1, 2), a slightly different event selection was used with respect to the T 2ν 1/2 analysis. If an event occurs with energy deposition in two detectors and the energy deposit in the detector where the decay took place is below the threshold for the anti-coincidence cut, the event contributes to the energy spectrum of the other detector. Therefore, when determining the total energy spectrum resulting from decays in one of the detectors, the energy spectra from all detectors in the array have to be taken into account. Such a selection has no impact on the detection efficiency for the Majoron process with n = 3 and 7 and 2νββ decay.  Fig. 2 Upper panel: experimental data (markers) and the best fit model (black histogram) for the golden data set. The contribution from 2νββ (green) and from the single background components are also shown. Lower panel: ratio between experimental data and the prediction of the best fit model. The green, yellow and red regions are the smallest intervals containing 68 %, 95 % and 99 % probability for the ratio assuming the best fit parameters, respectively [40].
array, for decays taking place in the active and dead part of detector α, becomes: with Φ α,0νχ AV,i,j (Φ α,0νχ DL,i,j ) giving the content of the i-th bin of the normalized energy distribution recorded with detector j for 0νββχ taking place in the active (dead) volume of detector α. Summing up the simulations of decays in all N det detectors results in the final model spectrum: For all four Majoron modes (n = 1, 2, 3, 7) only lower limits on the half-life can be given. They were obtained from the 90 % quantiles of the marginalized posterior distributions. These lower limits for T 0νχ 1/2 , not taking into account the systematic uncertainties, are in units of 10 23 yr: >4.4, >1.9, >0.9, and >0.4 for n = 1, 2, 3, and 7, respectively. The respective half-life of the 2νββ process derived from this analysis amounts to in units of 10 21 yr: 1.96±0.03 stat , 1.97±0.03 stat , 1.98±0.03 stat , and 1.99±0.03 stat . Within the uncertainties coming from the different background models and the different data sets of the two analyses, the derived T 2ν 1/2 values are in agreement (<1 σ) with that discussed in section 6.3.

Systematic Uncertainties
The systematic uncertainties were divided into the three categories (i) detector parameters and fit model, (ii) MC simulation, and (iii) data acquisition and selection.
(i) detector parameters and fit model Uncertainties from the fitting procedure were folded into the posterior distribution of T 0νχ 1/2 with a MC approach. Each source of uncertainty is described by a probability distribution. The fitting procedure was repeated 1000 times, each time drawing a random number for each source of uncertainty according to its probability distribution: -Material screening measurement results were used to constrain the minimum number of events expected from close and medium distant sources of the 214 Bi and 228 Th decays. Gaussian distributions describing these lower limits used in the fit were derived from the mean and standard deviations of the screening measurements. For details refer to Ref. [25]. -As for the T 2ν 1/2 analysis, the standard fit uses a bin width of 30 keV for the data and MC energy spectra. In order to determine the systematic uncertainty related to binning effects the bin width was sampled uniformly from 10 keV to 50 keV.
-Uncertainties on the active volume fractions enter the model in several ways. On the one hand, the MC energy spectra for all internal sources, that is for 2νββ, 0νββχ, 60 Co, and 68 Ga decays, are affected, as the fraction of decays taking place in the active and dead part of the detectors changes with changing f AV . On the other hand, the uncertainty on the active volume fraction also plays a role for the shape of the energy spectrum due to 42 K decays on the n + surface. Larger f AV means thinner n + dead layer and thus the possibility of an increased contribution from the electrons to the spectrum. For smaller f AV and thicker n + dead layer, their contributions are expected to be reduced. The active volume fraction for each detector was sampled from a Gaussian distribution with mean and standard deviation according to Table 1. For the coaxial detectors, the partial correlations of the uncertainty were taken into account. The simulated spectra of the internal sources as well as of the 42 K decays on the n + surface are composed according to the sampled active volume fractions. -The uncertainty on the fraction of enrichment in 76 Ge of the germanium that constitutes the detec-tors plays a role when converting the number of events attributed to 0νββχ into T 0νχ 1/2 . The probability distribution of f 76 for each detector is given by a Gaussian function with mean values and standard deviations as listed in Table 1.
-The data does not allow the resolution of the ambiguity regarding the exact positions of the near and medium distant sources. The 214 Bi decays serves as a representative in order to estimate the impact of this uncertainty. Their near position is represented by decays in the holders, in the mini-shroud or on the n + surface of the detectors, each having a probability of 1/3 in the sampling process. The medium distant position is represented by decays in the radonshroud or in the LAr, having a probability 1/2 in contrast. -Extensive studies of the characteristics of the BEGe diodes suggest the presence of a transition layer between the region where the detector is fully efficient and the external dead region [36,37]. An uncertainty as high as ±0.5 % on the lower limits of T 0νχ 1/2 is estimated for this effect in the case of the BEGe detectors. This uncertainty was folded into the total marginalized posterior distribution a posteriori. The corresponding uncertainty for the coaxial detectors is estimated to be negligible.
The marginalized posterior distributions for T 0νχ 1/2 derived from each of the 1000 individual fits were summed up. The resulting total marginalized posterior distribution accounts for the statistical as well as for the listed systematic uncertainties related to the fit model.
As for the T 2ν 1/2 analysis, the uncertainties on the active volume fractions and on the enrichment fractions are major contributions to the total uncertainty on the limits for T 0νχ 1/2 . However, the largest source of uncertainty is the composition of the fit model and the individual background contributions. In the case of n = 1, a fit with a bin width of 50 keV weakens the limit by ≈ 16 % compared to the standard fit, while the result for T 2ν 1/2 is not affected at all. The stability of the T 2ν 1/2 results shows the validity of the fit. The use of the alternative close and medium distant source positions for 214 Bi decays leads to maximal variations of +8.3 −12.6 % of the limit on T 0νχ 1/2 . (ii) MC simulation As in the case of the T 2ν 1/2 measurement, a total MC simulation uncertainty of 2.2 % has to be taken into account for effects related to the geometry implementation and particle tracking. It is folded into the total marginalized posterior distributions. No effect on the lower limits is observed for any of the spectral modes.
(iii) Data acquisition and selection Table 3 Experimental results for the limits on T 0νχ 1/2 of 76 Ge for the Majoron models given in Refs. [7,42,43,44]. The first section considers lepton number violating models (I) allowing 0νββ decay, while in the second section lepton number conserving models (II) are listed, where 0νββ decay is not allowed. The first column gives the model name, the second the spectral index, n, the third the information on whether one Majoron, χ, or two Majorons, χχ, is emitted, the fourth if the Majoron is a Goldstone boson, the fifth provides its lepton number, L, the sixth the experimental limit on T 0νχ 1/2 of 76 Ge obtained in this analysis. The nuclear matrix elements, M 0νχ , the phase space factor, G 0νχ , and the resulting effective coupling constants, g , are given in the seventh, eighth and ninth columns, respectively. The limits on T 0νχ 1/2 of 76 Ge for the Majoron models and g correspond to the 90 % quantiles of the marginalized posterior probability distribution. For the case of n = 1, the nuclear matrix element, M 0νχ , from Refs. [45,46,47,48,49,50,51] and the phase space factor, G 0νχ , from Ref. [52] are used for the calculation of g . The given range covers the variations of M 0νχ in these works. For n = 3 and 7, g is determined using the matrix elements and phase space factors from Ref. [42]. The results for 0νββχ (n = 3, 7) account for the uncertainty on M 0νχ . For n = 2, only the experimental upper limit is given. The uncertainty from data acquisition and selection is estimated to be below 0.1 % and does not alter the derived limits on T 0νχ 1/2 .

Results and Discussion
Fig . 3 shows the global model for the case of spectral index n = 1 together with the energy spectra for both the coaxial and the BEGe data sets. The contributions from the background contaminations, from the 2νββ decay only, and the combined spectra from the background contaminations and 2νββ decay are drawn separately. The 35868 events in the data spectrum of the golden data set were matched with 35834 events in the bestfit model for n = 1. Of those events, in the best fit, 54.5 are attributed to 0νββχ. For the BEGe data set, the best-fit model contains 5081.4 counts for the 5035 measured events. In this fit, 7.8 events are attributed to 0νββχ decay. The limit of T 0νχ 1/2 at 90 % C.I. derived from the fit is also drawn (green histogram). The upper limits at 90 % C.I. for the remaining three modes are reported for illustrative purpose (blue histogram for n = 2, orange for n = 3 and red for n = 7). The maximum of the corresponding distributions shifts to higher energy with the diminishing of the spectral index n. The resulting lower limits on T 0νχ 1/2 , determined as the 90 % quantiles of the posterior probability distributions and taking into account all uncertainties related to the fit model, are (in units of 10 23 yr): >4.2, >1.8, >0.8 and >0.3 for n = 1, 2, 3 and 7, respectively. The results are summarized in Table 3 for the different Majoron models.
The limits on T 0νχ 1/2 presented here are the most stringent limits obtained to date for 76 Ge. The limits for n = 1 and n = 3 are improved by more than a factor six [9], the limit for n = 7 is improved by a factor five [8] compared to previous measurements. The limit for the mode with n = 2 is reported here for the first time.
From the lower limits on T 0νχ 1/2 , upper limits on the effective neutrino-Majoron coupling constants g for the models with n = 1, 3 and 7 can be calculated using the following equations: and for single and double Majoron emission, respectively. The matrix element for the models with n = 1 (IB, IC and IIB) are taken from Refs. [45,46,47,48,49,50,51], whereas the phase space factor is that of Ref. [52]. The matrix elements for the models with n = 3 (ID, IE, IIC, IID, IIF) and with n = 7 (IIE) as well as the corresponding phase space factors are taken from Ref. [42]. The results for the upper limits on g are also Best-fit model and data energy spectrum for the coaxial and the BEGe data sets for the case of spectral index n = 1. The contributions from 2νββ decay and the background contributions are shown separately. The best-fit model does not contain the contributions from 0νββχ. The smallest interval of 68 % probability for the model expectation is indicated in grey. Also shown is the upper limit for 0νββχ decay with n = 1 as determined from the 90 % quantile of the marginalized posterior probability for 1/T 0νχ 1/2 . For illustrative purpose, also the upper limits at 90 % C.I. of the other three spectral indices n = 2, 3, 7 are reported.
shown in Table 3. The coupling constants allow a comparison with other isotopes. The best limits on 0νββχ decay of isotopes other than 76 Ge have been obtained for 100 Mo [10] and 136 Xe [14]. When comparing with the case of 100 Mo, it becomes obvious that the limits on T 0νχ 1/2 determined in the present analysis are about one order of magnitude more stringent, for the case of n = 7 even two orders of magnitude. However, due to the differences in the matrix elements and the phase space factors, the resulting limits on g from 100 Mo and 76 Ge are comparable. The limits for g derived from 136 Xe are a factor of two to five more stringent due to the higher limits that had been measured for T 0νχ 1/2 .

Conclusions
Phase I of the Gerda experiment, located at the INFN Laboratori Nazionali del Gran Sasso (LNGS) in Italy, has been executed between November 2011 and May 2013. Utilizing the collected exposure of Phase I, an improved result of the half-life of the 2νββ process in 76 Ge was obtained and new limits for the half-lives of the Majoron-emitting double beta decays were produced.
Thanks to the extremely low background level in the Gerda experiment, with a signal-to-background ratio of 3:1 in the 570-2039 keV interval and a refined background model, the measurement has an unprecedented precision (<5 %) with respect to previous experiments using 76 Ge. The new result is in good agreement with the one derived from a smaller data set with 5 kg·yr exposure [41]. The inclusion of more components into the reference background model results in a slight increase of the best estimate for T 2ν 1/2 . Majoron emission processes were searched for in the energy spectra using an exposure of 20.3 kg·yr. The analysis was performed for all four possibilities of the spectral index n (n = 1, 2, 3, and 7). No indication for a contribution of 0νββχ was found in any of the cases. Lower limits on the half-lives, T 0νχ 1/2 , were determined from the quantiles of 90 % probability of the marginalized posterior probability distributions. The results constitute the most stringent limits on T 0νχ 1/2 of 76 Ge obtained to date. For the standard mode (n = 1), the lower limit is determined to be: T 0νχ 1/2 > 4.2 · 10 23 yr.