AMoRE Collaboration

Experiments with solar and atmospheric neutrinos [1, 2] have found that neutrinos are massive. Various oscillation experiments [3, 4, 5, 6] have measured the three mixing angles and two mass differences. Although the absolute masses of neutrinos have not yet been measured, we know they are very small, less than 1 eV/$c^{2}$ based on measurements ofthe tritium $\beta$-spectrum endpoint [7, 8].

The small neutrino masses and the absence of right-handed neutrinos in the Standard Model motivated the introduction of Majorana masses for neutrinos, as opposed to charged leptons that have Dirac masses. Small Majorana neutrino mass can be generated via the seesaw mechanism in which the masses of active neutrinos are suppressed by heavy right-handed sterile neutrinos [9], and the mass terms violate lepton number conservation [10].

The currently well-established method for determining if neutrinos are Majorana fermions is to search for neutrinoless double beta $\rm(0\nu\beta\beta$) decay of nuclei [11]. Regardless of the underlying mechanism, the observation of $\rm 0\nu\beta\beta$ decay would prove that lepton number is violated [12], which is necessary for the seesaw mechanism and leptogenesis [13]. The amplitude of $0\nu\beta\beta$ decay is proportional to the effective Majorana mass, defined by the charged-current couplings of Majorana neutrinos. A $0\nu\beta\beta$ decay experiment requires monitoring a large mass of isotopes with a detector having an ultralow background and a high energy resolution to reach the true effective Majorana mass. In spite of more than seventy years of experimental efforts [14], no signal for $0\nu\beta\beta$ decay has been observed. The best half-life limit has been established by KamLAND-Zen for $\rm{}^{136}Xe$ with $T_{1/2}^{0\nu}>3.8\times 10^{26}$ yr, limiting the effective Majorana mass to $m_{\beta\beta}$$<$(28–122) meV [15]. For other isotopes, the $T_{1/2}^{0\nu}$ ($m_{\beta\beta}$) limits are: $\rm 3.3\times 10^{25}$ yr (90–305 meV) for $\rm{}^{130}Te$ by CUORE using cryogenic techniques [16], $\rm 1.8\times 10^{26}$ yr (79–180 meV) for $\rm{}^{76}Ge$ by GERDA using high purity germanium detectors [17], and $\rm 1.8\times 10^{24}$ yr (280–490 meV) for $\rm{}^{100}Mo$ by CUPID-Mo using cryogenic scintillating crystal detectors [18, 19].

AMoRE aims to search for the $0\nu\beta\beta$ decays using molybdate-based crystal detectors operating at milli-Kelvin (mK) temperatures [20, 21], similar to CUPID-Mo. AMoRE intends to achieve quasifree of background, defined as having less than one count in the region of interest (ROI) near $Q_{\beta\beta}$=3034.4 keV for the five-year experiment with about 100 kg of ¹⁰⁰Mo isotope.

The AMoRE project is progressing in three phases. The first two, AMoRE-Pilot and AMoRE-I, were completed at the 700-meter-deep Yangyang Underground Laboratory (Y2L). AMoRE-II will run at the newly built Yemilab for five years. We have reported a half-life limit of $T_{1/2}^{0\nu}=3.0\times 10^{23}$ yr at a 90% confidence level (CL) using six ^48deplCa¹⁰⁰MoO₄  crystals in the AMoRE-Pilot stage [22, 23, 24]. To confirm the performance and long-term stability of the detection system and determine the background level achievable with the existing setup at Y2L, we operated AMoRE-I from 2019 to 2023. This report describes the experimental setup, data analysis, and a new half-life limit of $\rm{}^{100}Mo$.

AMoRE-I was carried out in the same cryostat that was used for AMoRE-Pilot [22, 23, 24, 25], but with a larger number of detectors. We made a few modifications to the detector modules and shielding enhancements. The AMoRE-I system used a two-stage temperature control system to maintain the detector tower at a constant temperature [26]. The datasets used for this report were taken at a temperature of 12 mK.

The AMoRE-I detector array consists of eighteen crystals with a total mass of 6.2 kg, including 3.0 kg of ¹⁰⁰Mo, comprising five $\rm Li_{2}{}^{100}MoO_{4}$ (LMO) crystals, and thirteen ^48deplCa¹⁰⁰MoO₄ (CMO) crystals, six of which were inherited from AMoRE-Pilot. The molybdenum in the crystals was enriched in ¹⁰⁰Mo to 95.7$\pm$0.2%. Each detector module consisted of a crystal surrounded by a Vikuiti reflector film, an MMC sensor connected to the crystal, and a light detector made of an absorber with an MMC sensor. The mass-spring vibration damping system [27] was removed to accommodate the increased number of detector modules. The stainless steel screws in the detector modules were replaced with copper and brass for lower radioactive contamination. Additionally, Si wafers with SiO₂ antireflection coating were used for the light absorber instead of Ge wafers [28]. One flat surface of the LMO crystal, upon which the gold phonon collector was evaporated, was ground with 1500-grit SiC sandpaper to enhance the gold film bond at the interface. Each crystal had a stabilization heater attached to the flat surface to measure and correct the gain drift of the heat signal [29].Pulses of fixed current and duration were injected through the heater every 10 seconds. Detailed descriptions of the detector module can be found in  [30, 26].

An additional 5 cm of lead was added outside of the cryostat vacuum chamber to further reduce high-energy environmental $\gamma$’s. Initially, aluminum plates were used to support the lead bricks, but these plates were found to be highly contaminated by $\rm{}^{228}Ra$ and were removed in May 2021. Ten additional muon counters made of plastic-scintillator panels were installed to extend the solid-angle coverage for the detector array.

Heat and light signals from the SQUID electronics were continuously digitized without an event trigger and stored using flash analog-to-digital converter modules with an 18-bit resolution for a 10-volt peak-to-peak dynamic range at a 100 kHz sampling frequency. Data acquisition (DAQ) with the full eighteen-crystal detector array started in December 2020 and ended in April 2023, with a 93% live time. About 78% of the DAQ live time was dedicated to physics measurement, while calibration and other commissioning data took up the remaining time. Thorium-containing welding rods were placed between the outer vacuum chamber and the external lead shielding to calibrate the energy scale every 4–8 weeks, depending on the detector’s stability.

To suppress noise, the heat signal of each detector was processed using a Butterworth bandpass filter for the offline analysis. For most detectors, the trigger thresholds were set below 100 keV. Two modules that suffered from a significant vibration or electric noise had higher thresholds. The amplitudes of the heat and light signals were determined by fitting template waveforms to the data using a least-squares fit [31]. The fit was applied to the filtered waveforms and the bandwidths of the Butterworth filters were optimized to provide the highest energy resolutions. A template signal was constructed by averaging 2615 keV $\gamma$ events accumulated in calibration runs. The rise and fall times were derived from the raw signals, defined as the time elapsed between specific points on the pulse. Unlike the AMoRE-pilot analysis, which used the time difference between the 10% and 90% levels of the signal, this analysis defined rise time (RT) variably for each detector, optimizing the pulse shape discrimination (PSD) power. One LMO detector had a noisy light signal, making it impossible to calculate the light signal amplitude accurately. Data from this detector were excluded from the analysis because of the limited PSD power. Examples of the RT versus signal amplitudes for a CMO and an LMO detector are shown in Fig 1.

The detector response change over time, mainly caused by temperature variations, influenced the pulse amplitudes and shapes of both physics and heater signals. Correlations between the heater signal’s RT and the amplitude of 2615 keV $\gamma$ events were determined with calibration data and used for drift correction in the corresponding dataset. Each dataset consists of a calibration run followed by physics runs before the next calibration.

The $\beta/\gamma$ energy scale was initially calibrated using four prominent $\gamma$ peaks at 511.0 keV, 583.2 keV, 911.2 keV, and 2614.5 keV in the calibration-run spectrum for each dataset. As in the previous analysis, signal amplitudes for each detector were all well described as a quadratic function of true energy values, without a constant term [32, 22]. Secondary energy calibration was performed using another quadratic function with a constant term, including more $\gamma$ peaks from the calibration data.

Owing to nonuniform, position-dependent detector responses, peaks in the energy spectra were asymmetric, showing tails on both the lower and upper sides. Among various peak-fitting functions, the Bukin function [33] was found to be well suited for the peaks in the calibration spectrum. The spectra around the 2615 keV peak in the calibration data for a CMO and an LMO detector are shown in Fig. 2. Shape parameters that were determined from fitting below 2615 keV were extrapolated to $Q_{\beta\beta}$ to estimate the $0\nu\beta\beta$ signal shape. The squares of energy resolutions ($h^{2}$) were well fitted to a quadratic function of energy. The energy resolution at $Q_{\beta\beta}$ differs for each detector, varying from 10 to 28 keV of full width at half maximum with a weighted average 13 keV. Other shape parameters for asymmetry ($\xi$), left tail ($\rho_{L}$), and right tail ($\rho_{R}$) varied slowly with energy, and their uncertainties were extrapolated using a linear function or were left constant, depending on each detector’s characteristics.

The $0\nu\beta\beta$ event-selection criteria for background reduction were divided into two categories. The first category was particle identification to remove the continuous $\alpha$ background around the ROI due to radioactive contamination on the crystal surface or surrounding materials. Two parameters were adopted for each detector: PSD using RT of the heat signal, and the light-to-heat ratio (L/H), which leverages the differences in scintillation quenching for $\alpha$ and $\beta/\gamma$ particles. Fig. 1 shows the RT and L/H versus energy for the calibration runs, background runs, and events that were coincident with muons for one CMO and one LMO detector. Here, the $\alpha$ events are clustered at higher energies and at smaller RT and L/H values, while most of the $\beta/\gamma$ events without muon coincidence lie below 2.6 MeV.

The median-absolute-deviation (MAD) values were determined using source data at the energy range of 465–2665 keV and fitted using a function of energy. Events within the 3 MAD bands for RT and L/H, denoted as the dashed-orange curves in Fig. 1 were selected as $\beta/\gamma$ events. The selection efficiencies at the calibration $\gamma$-peak energies were determined by comparing the event counts at each peak before and after the selection using the Bukin function on top of a simple exponential or linear background to extract the peak from the resulting spectra. These $\beta/\gamma$ selection efficiencies, when extrapolated to $Q_{\beta\beta}$, varied among detectors and ranged between 80% and 95% with uncertainties on the 1% level. Generally, CMO detectors showed a much better discrimination power than LMO, for both RT and L/H parameters, due to the higher light yield and the correlation of light yield with the pulse shape [34, 35]. The alpha backgrounds for the six CMO detectors used in AMoRE-Pilot were reported in [23]. The activities of $\rm{}^{226}Ra$ ($\rm{}^{228}Th$) for the AMoRE-I dataset were determined by analysis of $\alpha-\alpha$ coincident events from the sequential decay of ²²²Rn and ²¹⁸Po (²²⁴Ra and ²²⁰Rn). These were 5–59 $\mu$Bq/kg (1–13 $\mu$Bq/kg) in the CMO detectors, except for one highly contaminated CMO crystal, and 1–2 $\mu$Bq/kg for both $\rm{}^{226}Ra$ and $\rm{}^{228}Th$ in the LMO detectors. Detailed background modeling of AMoRE-I is in progress.

The second event-selection category involved a series of anticoincidence cuts. First, events occurring within 10 ms after a muon candidate event were rejected. About ten thousand muon candidate events were registered daily, resulting in an inefficiency of only 0.1% with negligible uncertainty. Secondly, we imposed a single-hit cut in which the signal candidate is the only hit that occurred within a 2 ms time window. Since the trigger rate of the background runs was typically $\sim$1 Hz, the efficiency of this single-hit requirement was about 99.8%. The final anticoincidence condition, called $\alpha$ tagging, targeted one of the major backgrounds in the ROI from the $\beta$ decay of ²⁰⁸Tl with $Q\sim$4.998 MeV. These background events follow an $\alpha$ event emitted in ²¹²Bi decay to ²⁰⁸Tl. Since the half-life of ²⁰⁸Tl is 3.053 minutes and the energy released in the $\alpha$ decay of ²¹²Bi is 6.207 MeV, all events within 20 minutes after an $\alpha$ event in the same crystal with 6.2 $\pm$ 0.1 MeV energy were rejected. Two detectors with exceptionally high $\alpha$ rates had a narrower $\alpha$ energy window of 40 keV [23]. The efficiency of this $\alpha$ veto was about 98% on average: 78% for the most $\alpha$-contaminated CMO and higher than 99% for LMOs. The cut efficiencies are summarized in Table 1.

Figure 3 shows the resulting background energy spectra accumulated for 8.02 kg$\cdot$year (or 3.89 kg${}_{\mathrm{{}^{100}Mo}}\cdot$year) exposure, following the step-by-step selections, overlaid with the calibration spectrum. The first dataset taken with the aluminum support plate, shown as the dark blue shaded area, represents about 14% of the total exposure but contains about 20% of the events in the energy range below 2.7 MeV, including about 60% of 2615 keV $\gamma$ events. Despite this, the background due to the contaminated aluminum plate is mainly confined to energies below 2.7 MeV, so these data were included in the $0\nu\beta\beta$ analysis. As shown in the previous studies [22], the relationship between the signal amplitude and energy differs for $\alpha$ events compared to $\beta/\gamma$ events due to pulse shape differences, which vary across detectors. Consequently, $\alpha$ energies, which were determined separately, are not reported in this work. As a result, the energy spectrum for all detectors and datasets shown in Fig. 3 has many $\alpha$ peaks above 4 MeV electron-equivalent energy that are distributed incoherently when they are calibrated with the functions for $\beta/\gamma$ events.

Events in the 2.7–3.6 MeV energy range were selected for the $0\nu\beta\beta$ study. The background energy distribution was approximated as a linear combination of flat and exponential background components. The exponential component includes high-energy $2\nu\beta\beta$-decay events and $\beta$ background from internal and surface contamination of the crystal [23]. The flat component is introduced to describe high energy $\gamma$s from neutron capture, external sources such as rock and radon in the air, and residual $\beta$ and $\alpha$ events due to incomplete rejection. A more detailed background modeling based on the measurements of detector components’ radioactivities is in progress.

Considering the background ($b$) and the $0\nu\beta\beta$ decay signal ($s$) can be expressed in terms of the decay rate ($\Gamma^{0\nu}$=$\ln 2/T^{0\nu}_{1/2}$), the number of ¹⁰⁰Mo nuclei ($N_{100}$), the detection efficiency ($\varepsilon$) for the $0\nu\beta\beta$ decay of $\rm{}^{100}Mo$ (shown in Table 1), and the DAQ live time ($t$) as $s$=$\varepsilon\Gamma^{0\nu}N_{100}t$. An unbinned likelihood function was constructed as follows:

where $i$ and $j$ are indices of the detector and its events, respectively, and $n_{i}$ is the number of observed events in the $i$th detector. The expected spectral shape was described by the probability density function, $f^{s+b}$, that includes a background model, with flat and exponential components, plus a signal peak in the form of a Bukin function for fully contained $0\nu\beta\beta$ events. Efficiency ($\varepsilon_{i}$) and peak shape parameters ($\bm{\nu}_{i}$) were treated as nuisance parameters with Gaussian priors ($\pi$).

The minimized negative log likelihood was profiled over the sensitive range of $\Gamma^{0\nu}$ from 0 to $\sim 10^{-23}$ year^-1, and the best fit was found at $\Gamma^{0\nu}=0.0$ year^-1, meaning that no event excess was found over the assumed background shape, as shown in Fig. 4. The corresponding limit on the half-life of the $0\nu\beta\beta$ decay of ¹⁰⁰Mo was evaluated at a 90% CL to be:

which extends the CUPID-Mo limit [19]. The limit matched well to the sensitivity determined by numbers of pseudodatasets with the given exposure, peak shapes and sideband background levels. The median $\pm$34% sensitivity of 90% CL limit is 2.8${}^{+1.0}_{-0.9}\times 10^{24}$ yr, as shown in Fig. 4.

The posterior analysis leads to a total background level of $b=0.025\pm 0.002$ counts/keV/kg/yr. The LMO detectors showed a slightly lower background rate of $0.021\pm 0.005$ counts/keV/kg/yr on average, while the background counting rate of the CMO detectors is $0.026\pm 0.003$ counts/keV/kg/yr. The total background rate in AMoRE-I was reduced by $\sim$15 times compared to that of AMoRE-pilot.

We calculated the effective Majorana mass within the theoretical framework of the light neutrino exchange model, incorporating the phase space factor [39, 40, 41] and nuclear matrix elements (NMEs) [42, 43, 44, 45, 46, 47, 48, 49]. The range of inferred upper limits on the effective Majorana mass is $m_{\beta\beta}$$<$(210–350) meV assuming an axial vector coupling constant $g_{A}=1.27$. If we include the first shell model calculation recently published for $\rm{}^{100}Mo$ [50], the range of limits extends to $m_{\beta\beta}$$<$ (210–610) meV. The lower limit is derived from an energy density functional considering the nuclear deformation and pairing fluctuation [44], while the upper limit is based on a shell model that explicitly includes the short-range correlations [50]. Upper limits on the effective Majorana mass from existing experimental data and NMEs are shown in Fig. 5.

AMoRE-II is under preparation at Yemilab [53] with a muon rate about a quarter of that at Y2L [54]. We have developed an LMO detector module with improved energy resolution and alpha background rejection [55, 56]. The background level for AMoRE-II is projected to be less than $10^{-4}$ counts/keV/kg/yr based on radioassay data and GEANT4 simulations [57]. The discovery sensitivity is projected to be approximately $\rm 4.5\times 10^{26}$ yr for five years of data collection.

Acknowledgements–This research is supported by Grants No. IBS-R016-D1 and No. IBS-R016-A2. It is also supported by the National Research Foundation of Korea (NRF-2021R1I1A3041453, NRF-2021R1I1A6A1A03043957) and the National Research Facilities & Equipment Center (NFEC) of Korea (No. 2019R1A6C1010027). We appreciate the support by the Ministry of Science and Higher Education of the Russian Federation (N121031700314-5), the MEPhI Program Priority 2030. The group from the Institute for Nuclear Research (Kyiv, Ukraine) was supported in part by the National Research Foundation of Ukraine under Grant No. 2023.03/0213. We thank the Korea Hydro and Nuclear Power (KHNP) Company for providing underground laboratory space at Yangyang and the IBS Research Solution Center (RSC) for providing high-performance computing resources. These acknowledgements are not to be interpreted as an endorsement of any statement made by any of our institutes, funding agencies, governments, or their representatives.