Observations of the neutrino oscillations suggest that the neutrinos are massive, which calls for extension of the Standard Model of particles and fields (SM). However, oscillation experiments cannot determine the neutrino mass and the neutrino mass hierarchy. One of the most promising tools to determine the absolute neutrino mass scale and the neutrino mass hierarchy, the nature of the neutrino (Dirac or Majorana particle?), to check the lepton number conservation is double beta (2β) decay of atomic nuclei, a process in which two electrons (or positrons) are emitted simultaneously and nuclear charge changes by two units: (A,Z)→(A,Z±2) [1] [2] [3] . The neutrinoless mode of the decay (0ν2β) violates the lepton number conservation law and it is possible if the neutrinos are Majorana particles (particle is equal to its antiparticle). Being a process beyond the SM, the 0ν2β decay has the potential to test the SM [4] [5] [6] . Moreover, the Majorana nature of the neutrino might shed light on the Universe baryon asymmetry problem [7, 8] .
The two-neutrino 2β decay (2ν2β) is a radioactive process allowed in the SM with the longest half-lives ever observed: 10 18 -10 24 yr. The 2ν2β − decay mode has been detected in several nuclides [9] . the STELLA (SubTErranean Low Level Assay) facility [33] at the LNGS. The first stage of the experiment with the 106 CdWO 4 detector gave the half-life limits on 2β processes in 106 Cd at level of ∼ 10 20 yr [21] . In the second stage the 106 CdWO 4 scintillator was installed between four HPGe detectors (with volume 225 cm 3 each) of the GeMulti HPGe γ spectrometer to detect γ quanta expected in the most of the 106 Cd decay channels, including the annihilation γ's emitted in decay modes with positron(s) emission (a simplified decay scheme of 106 Cd is presented in Fig. 1 ). The experiment improved the 106 Cd half-life limits to the level of T 1/2 ≥ (10 20 − 10 21 ) yr [34]. In the third stage, described in the present report, the 106 CdWO 4 detector was running in coincidence (anti-coincidence) with two large volume CdWO 4 crystal scintillators in a close geometry to increase the detection efficiency to γ quanta expected to be emitted from the 106 CdWO 4 crystal in the double beta decay processes in 106 Cd. Preliminary results of the experiment stage were reported in [35] .
The 106 CdWO 4 crystal scintillator of roughly cylindrical shape (approximate sizes 27 mm ×50 mm, mass 215.4 g) was viewed by a 3 inches low radioactive photo-multiplier tube (PMT) Hamamatsu R6233MOD through a lead tungstate (PbWO 4 ) crystal light-guide ( 40 mm ×83 mm). The PbWO 4 crystal has been developed from the highly purified [37] archaeological lead [38] . Two CdWO 4 crystal scintillators 70 mm ×38 mm include a cylindrical cut-out to house the 106 CdWO 4 crystal. They were viewed by two 3 inches low radioactive PMTs EMI9265B53/FL through light-guides glued in two parts: low radioactive quartz ( 66 mm ×100 mm, close to the CdWO 4 scintillators) and optical quality polystyrene ( 66 mm ×100 mm). A schematic of the set-up is shown in Fig. 2 . The detector system was surrounded by four high purity copper bricks (referred hereinafter as "internal copper") and by layers of high purity copper (1) is viewed through PbWO 4 light-guide (2) by photo-multiplier tube (3) . Two CdWO 4 crystal scintillators (4) are viewed through light-guides glued from quartz (5) and polystyrene (6) by photo-multiplier tubes (7). The detector system was surrounded by passive shield made from copper, lead, polyethylene and cadmium (not shown). Only part of the copper details (8, "internal copper"), used to reduce the direct hits of the detectors by γ quanta from the PMTs, are shown.
(11 cm, referred hereinafter as "external copper"), low radioactive lead (10 cm), cadmium (2 mm) and polyethylene (10 cm) to reduce the external background. The inner volume of the set-up with the detector system was continuously flushed by high-purity nitrogen gas to remove environmental radon. The grade of the high-purity N 2 gas is at least 5.5 for what concerns the presence of other possible gases. However, the possible presence in trace of Radon gas in the Nitrogen atmosphere inside the copper box, housing the detector, has been checked with another set-up, by searching for the double coincidences of the γ-rays (609 and 1120 keV) from 214 Bi Radon daughter. The obtained upper limit on the possible Radon concentration in the high-purity Nitrogen atmosphere has been measured to be: < 5.8 × 10 −2 Bq/m 3 (90% C.L.) [39] . Photographs of the detector system are shown in Fig. 3 .
An event-by-event data acquisition system based on a 100 MS/s 14 bit transient digitizer (DT5724 by CAEN) recorded the amplitude, the arrival time and the pulse shape of each event. To reduce the data volume due to presence in the 106 CdWO 4 crystal of 113 Cd and 113m Cd β active nuclides [21, 32] , the energy threshold for the set-up was set at level of ≈ 510 keV for the anti-coincidence mode, while the energy threshold of the 106 CdWO 4 detector in the coincidence with the CdWO 4 counters was ≈ 200 keV. The energy thresholds of the CdWO 4 counters were ≈ 70 keV. The energy scale and the energy resolution of the detectors were measured with 22 Na, 60 Co, 133 Ba, 137 Cs, and 228 Th γ sources at the beginning, in the middle, and at the end of the experiment.
The energy resolution of the 106 CdWO 4 detector for the total exposure can be described by the function FWHM = 6.85 × E γ , where FWHM (full width at half maximum) and E γ are given in keV. The poor energy resolution of the enriched detector (despite excellent optical properties of the material [32] ) is caused by the elongated shape of the enriched scintillator that results in a rather low and non-uniform light collection, and by using of not perfectly transparent PbWO 4 crystal light-guide. The performance of the CdWO 4 counters is substantially better. The energy spectra accumulated by one of the counters with 22 Na, 60 Co and 228 Th γ sources are presented in Fig. 4 . The energy resolution of the counters was estimated by using the results of the three energy calibration campaigns as FWHM = a × E γ with the (3), "internal copper" bricks (4), "external copper" bricks (5), lead bricks (6), polyethylene shield (7) . The copper, lead and polyethylene shields are not completed. coefficient a equal to 2.97 and 3.13 for the two detectors. The resolution formulas take into account also energy scale shifts during the data taking over the experiment.
Energy spectra of 22 Na source were simulated by the EGSnrc code [40] . The data measured with 22 Na source without coincidence selection and in coincidence with energy 511 keV in at least one of the CdWO 4 counters is compared with the simulated distribution in Fig. 5 . The experimental data is in a reasonable agreement with the results of simulations.
A distribution of the 106 CdWO 4 detector pulses start positions relative to the CdWO 4 signals with energy 511 keV is shown in Inset of Fig. 5 . The time resolution of the detector system is rather high (the standard deviation of the distribution is 16 ns) due to the fast rise time of the CdWO 4 scintillation pulses.
The difference in CdWO 4 scintillation pulse shape for β particles (γ quanta) and α particles can be used to suppress the background caused by α radioactive contamination of the detector due to the residual contamination in 232 Th and 238 U with their daughters. The mean time method was applied to the data to discriminate signals of different origin by pulse shape. For each signal f (t), the numerical characteristic of its shape (mean time, ζ) was defined by using the following equation:
where the sum is over the time channels k, starting from the origin of signal up to 35 µs; f (t k ) is the digitized amplitude (at the time t k ) of a given signal. The energy dependence of the parameter ζ and its standard deviation (the distributions of ζ for β particles (γ quanta) and α particles are well described by a Gaussian function) was determined by using the data of the calibration measurements with 228 Th gamma source. The obtained parameters were then used to discriminate β (γ) events from α events in the data of the low-background experiment. We refer reader to our previous works [21, 34] where the pulse-shape discrimination (PSD) method was described in detail. By using the PSD the α events were statistically separated from γ(β) events. In addition the method discarded from the data events of the 212 Bi -212 Po sub-chain from the 232 Th family (due to the short decay time of 212 Po ≈ 0.3 µs these decays are treated by the data acquisition system as a single event), PMT noise, pile-ups of signals in the 106 CdWO 4 detector, 106 CdWO 4 plus PbWO 4 events, etc. The results of the PSD method application to the background data gathered for 26033 h in the low-background set-up is shown in Fig. 6 . The mean time method reduced the background mainly in the energy region (800-1300) keV (by a factor ∼ 1.6) where α events of the 232 Th and 238 U with their daughters are expected.
Further reduction of the background counting rate (by a factor ∼ 1.3 in the energy interval (1000-3000) keV) was achieved by exploiting the anti-coincidence with the CdWO 4 counters. The background was suppressed significantly by selection of events in the 106 CdWO 4 detector in coincidence with the event(s) in at least one of the CdWO 4 counters with the energy release E = 511 ± 2σ keV (by a factor ∼ 17 in the same energy interval; here σ is the energy resolution of the CdWO 4 counters for 511 keV γ quanta), and by selection of events in coincidence with the events in both the CdWO 4 counters with the energy E = 511 ± 2σ keV (by a further factor ∼ 42). The stages of the background spectra reduction are presented in Fig. 6 .
The counting rate of the 106 CdWO 4 detector below the energy of ∼ 0.8 MeV is mainly caused by the β decay of 113 Cd with the energy Q β = 323.83 (27) keV [22] and of 113m Cd (Q β = 587.37 (27) keV [22, 41] ). A background model to describe the experimental data after the 113m Cd β spectrum was constructed from distributions of "internal" (radioactive contamination of the 106 CdWO 4 crystal) and "external" (radioactive contamination of the set-up details) sources. The equilibrium of the 238 U and 232 Th chains in all the Table 1 . Radioactive contamination (mBq/kg) of the materials of the low-background set-up estimated by using the fit of the energy spectra presented in Fig. 7 . Upper limits are given at 68% C.L. 
materials is assumed to be broken 1 [47] .
The following "external" sources were simulated in the materials of the set-up: 56 Co and 60 Co in the internal copper bricks.
The background components were simulated using the EGSnrc package with initial kinematics given by the DECAY0 event generator [48] . The distribution of residual α particles of 232 Th and 238 U with their daughters was constructed from the experimental data by using the pulse-shape analysis.
The simulated models were used to fit the energy spectra of γ and β events in anti-coincidence with the CdWO 4 counters and in coincidence with event(s) in at least one of the CdWO 4 counters with the energy release E = 511 ± 2σ keV. The data were fitted in the energy intervals (940-4000) keV (anti-coincidence data) and (240-3940) keV (coincidence with 511 keV). The fit quality is reasonable (χ 2 = 457 for 235 degrees of freedom). The results of the fit and the main components of the background are shown in Fig. 7 .
The fit allowed to estimate limits on radioactive contamination of the materials of the low-background set-up. The data are presented in Table 1 .
There are no peculiarities in the experimental data that could be ascribed to 2β processes in 106 Cd. Lower limits on the half-life of 106 Cd relatively to different 2β decay channels and modes can be estimated by using the following formula: 1 Secular equilibrium in the 232 Th and 238 U decay families (when activities of daughter nuclides are equal to the activity of their parent nuclide) is typically broken in almost all materials due to physical or chemical processes utilized in the material production (see, e.g., [42] [43] [44] . The main components of the background are shown: the distributions of internal contaminations ("int 40 K", "int 232 Th", and "int 238 U") and external γ quanta ("ext γ"), residual α particles in the 106 CdWO 4 crystal (α), cosmogenic 56 Co and 60 Co in the copper shield details, the 2ν2β decay of 116 Cd. The excluded distributions of the 0ν2EC decay of 106 Cd to the ground state of 106 Pd with the half-life T 1/2 = 6.8 × 10 20 yr are shown by red solid line.
lim
where N is the number of 106 Cd nuclei in the 106 CdWO 4 crystal (2.42 × 10 23 ), η det is the detection efficiency for the process of decay (calculated as a ratio of the events number in a simulated distribution to the number of generated events), η sel is the selection cuts efficiency (selection by PSD, time coincidence, energy interval), t is the time of measurements, and lim S is the number of events of the effect searched for, which can be excluded at a given confidence level (C.L.). The responses of the detector system to different modes and channels of 106 Cd double beta decay were simulated using the EGSnrc package with initial kinematics given by the DECAY0 event generator. About 5 × 10 6 events were generated for each decay channel. Different data were analyzed to estimate limits on the 2β processes in 106 Cd. Fit of the anti-coincidence spectrum by the above described model plus a simulated distribution of the 0ν2EC decay of 106 Cd to the ground state of 106 Pd returns the area of the distribution (205 ± 99) counts that is no evidence for the effect searched for. According to [49] we took 367 events as lim S at 90% C.L. 2 The detection efficiency for the decay was simulated as η det = 0.522. Taking into account the selection cut efficiency due to application of the PSD to select γ and β events η sel = 0.955, we got a lower limit on the half-life of 106 Cd relative to the 0ν2EC decay to the ground state of 106 Pd T 1/2 ≥ 6.8 × 10 20 yr (the excluded distribution of the 0ν2EC decay is shown in Fig. 7) . The limit is slightly worse than the one obtained in the previous stage of the experiment (T 1/2 ≥ 1.0 × 10 21 yr [21] , see also Table 2 ).
Fit of the 106 CdWO 4 detector data in coincidence with signal(s) in the CdWO 4 counters by the above described background model was more sensitive to the most of the modes and channels of the decay searched for. An example of such an analysis for the 0νECβ + and 0ν2β + decays of 106 Cd to the ground state of 106 Pd by using the data measured with the 106 CdWO 4 detector in coincidence with 511 keV events in at least one of the CdWO 4 counters is shown in Fig. 8 . The selection cuts efficiency, e.g., for the 0νECβ + process was calculated to be η sel = 0.909 as a product of the PSD to select γ and β events in the interval ±2σ of the mean time values (0.9546), the time coincidence efficiency in the interval ±3σ (0.9973), the energy interval ±2σ to select 511 keV events in the CdWO 4 counters (0.9545). The data on the efficiencies, values of lim S and the obtained half-life limits are given in Table 2 .
Another example is search for 0ν2EC transition of 106 Cd to the 2718 keV excited level of 106 Pd (considered as one of the most promising decay channels from the point of view of a possible resonant process [14] ). The search was realized by analysis of the 106 CdWO 4 detector data in coincidence with event(s) in at least one of the CdWO 4 counters in the energy interval (1046 − 1.5σ)-(1160 + 1.7σ) keV. The interval should contain two intensive γ quanta with energies 1046 keV and 1160 keV expected in the decay searched for (see the decay scheme in Fig. 1) . The spectrum and its fit, consisting of the background model and excluded distribution of the resonant process searched for, is presented in Fig. 9 .
The highest sensitivity to several decay channels with positron(s) emission was achieved by using the data gathered by the 106 CdWO 4 detector in coincidence with 511 keV annihilation γ quanta in both of the CdWO 4 counters thanks to a rather high detection efficiency of the CdWO 4 counters and a very low background counting rate (see Fig. 10 ). However, the fit of the spectrum by the background components is not reliable enough due to a very low statistics of the data. Thus, the method of comparison of the measured background with the expected one was applied for the analysis. The expected background was estimated from the results of the fit shown in Fig. 7 . There are 54 counts in the whole spectrum, while the estimated background is 55.3 counts confirming a correct background modelling. In the energy interval 2 In the present work all the limits are given with 90% C.L. Only statistical errors coming from the data fluctuations were taken into account in the estimations of the lim S values, and systematic contributions have not been included in the half-life limit values. (250-1000) keV the measured background is 33 counts, while the estimated one is 37.4 counts that leads to lim S = 6.7 counts in accordance with the recommendations [49] . Taking into account the detection and the selections efficiencies for the 2νECβ + decay of 106 Cd to the ground state of 106 Pd (0.040 and 0.703, respectively) one can get a half-life limit T 1/2 = 2.1 × 10 21 yr that is about two times higher than the limit (T 1/2 = 1.1 × 10 21 yr) obtained in the previous stage of the experiment [34] . Limits on other 2β decay processes in 106 Cd were obtained in a similar way. They are presented in Table 2 , where the results of the most sensitive previous experiments are given for comparison.
A limit on effective nuclear matrix elements for the 2νECβ + decay of 106 Cd to the ground state of 106 Pd could be estimated by using the calculations of the phase-space factors for the decay [50, 51] with the formula 1/T 1/2 = G 2νECβ + × |M e f f | 2 . The effective matrix nuclear element M e f f is expressed by M e f f = g 2 A × M 2νECβ + , where g A is the axial-vector coupling constant, M 2νECβ + is nuclear matrix element. An upper limit on the value of the effective matrix nuclear element for the process can be estimated as M e f f ≤ (0.80 − 0.82).
The half-life limit on the 2νECβ + decay of 106 Cd to the ground state of 106 Pd, T 1/2 ≥ 2.1 × 10 21 yr, approaches the region of the theoretical predictions that are in the range 10 21 − 10 22 yr [10, [52] [53] [54] [55] . The sensitivity to the double beta decay processes in 106 Cd is expected to be improved in the currently running experiment with reduced background thanks to utilization of ultra-radiopure PMTs, longer quartz light-guides for the CdWO 4 counters, a more powerful passive shield of the detector system. Also the energy resolution of the 106 CdWO 4 detector was improved, roughly by a factor ∼ 1.8, thanks to replacement of the PbWO 4 light-guide by a plastic scintillator light-guide with a substantially better optical transmittance. This replacement became possible due to an extremely low radioactive contamination of the specially developed R11065-20 MOD Hamamatsu PMT [56] used for the 106 CdWO 4 detector. 
The experiment to search for double beta decay of 106 Cd with enriched 106 CdWO 4 scintillator in coincidence with two large volume CdWO 4 scintillation counters was performed at the Gran Sasso underground laboratory of INFN (Italy). New improved limits are set on the different channels of 106 Cd double beta decay at level of 10 20 − 10 22 yr. The new improved limit on half-life of 106 Cd relative to the 2νECβ + decay was estimated as T 1/2 ≥ 2.1 × 10 21 yr. The sensitivity is within the region of the theoretical predictions for the decay probability that are in the range of T 1/2 ∼ 10 21 − 10 22 yr. A new improved limit was set for the resonant neutrinoless double-electron capture to the 2718 keV excited level of 106 Pd as T 0ν2K 1/2 ≥ 2.9 × 10 21 yr. The next stage of experiment is running at LNGS in the DAMA/R&D set-up with an improved sensitivity to all the decay channels thanks to reduction of the background approximately by a factor 3-5 with utilization of ultra-radiopure PMTs, longer quartz light-guides for the CdWO 4 counters, a more powerful passive shield of the detector system. The energy resolution of the 106 CdWO 4 detector was improved too thanks to replacement of the PbWO 4 light-guide by a plastic scintillator light-guide with a substantially better optical transmittance. As a result, the sensitivity to the 2νECβ + decay of 106 Cd is expected to be high enough to detect the process with the half-life at level of ∼ (0.5 − 1) × 10 22 yr over 5 yr of measurements. 