Simulation Tool Development and Sensitivity Analysis of ¹⁶⁰Gd Double Beta Decay Search by the PIKACHU Project

We considered two types of background sources: nuclides in ²³²Th, ²³⁸U, and ²³⁵U decay chains contaminated in GAGG, and ⁴⁰K present in both GAGG and photomultiplier tube (PMT).

It has been confirmed through background measurement in Kamioka underground facility that high-purity GAGG contains ²³²Th, ²³⁸U, and ²³⁵U decay chains PIKACHU . On the other hand, the investigations using a high purity germanium detector demonstrated that the purity of crystal’s raw materials is higher than that of high-purity GAGG. This suggests that radioactive impurities may also be introduced from the environment during the crystal growth process.
In GEANT4, we generated the stationary parent nuclide for each decay chain, corresponding to ²³²Th, ²³⁸U, and ²³⁵U, uniformly within the volume of the crystal model. Each parent nuclide was generated 1,500,000 times in this study. GEANT4 then simulated their chain decays until the parent nuclide reached the most stable nuclide in each chain and the behavior of resulting radiations in the crystal model. The energy losses of the radiations in the crystal model were recorded, along with the names of the decayed nuclide and the daughter nuclide. However, the recorded energy losses do not reflect the effect of quenching or the detector’s energy resolution. These effects are artificially incorporated into analysis as follows.
Quenching is the effect in which $\alpha$-rays are detected with lower energies than their original energies. The ratio of the $\alpha$-ray energy calibrated by $\gamma$-rays to the original $\alpha$-ray energy (the $\alpha$/$\gamma$ ratio) is known to depend proportionally on the original $\alpha$-ray energy in GAGG quenching . Therefore, we approximated the correlation between the $\alpha$-ray energy before and after quenching using a quadratic function, as shown in Fig. 1 (b). The red line represents the quadratic fit function ($y=\mathrm{p_{0}}+\mathrm{p_{1}}\,x+\mathrm{p_{2}}\,x^{2}$) to the blue filled points with error bars. The blue points are obtained by fitting the measured $\alpha$-ray spectrum of high-purity GAGG PIKACHU with Gaussian functions (red lines) in Fig. 1 (a). Based on the curve in Fig. 1 (b), we corrected the energy losses of the $\alpha$-rays calculated by GEANT4 to the $\gamma$-equivalent energies. Then, by applying the detector’s resolution to the corrected energy losses, we reproduced the measured $\alpha$-ray spectrum. The resolution was measured using $\gamma$-ray sources and approximated with an offset of 0.148 MeV in Ref. PIKACHU . For this reason, the resolution cannot be applied to $\gamma$-equivalent energies less than 0.148 MeV in the model spectra presented in this paper. However, the effect of this offset is negligible because it is sufficiently lower than both the threshold of the data acquisition and the energy region of interest in our analysis.
Radioactive nuclides are present in GAGG in a state of radioactive equilibrium in relation to half-lives as shown in Table 1. We classified the ²³⁸U decay chain into five sections: ²³⁸U_up, ²³⁴U, ²³⁰Th, ²³⁸U_mid, and ²³⁸U_down, with ²³⁸U, ²³⁴U, ²³⁰Th, ²²⁶Ra, and ²¹⁰Pb as the parent nuclides, respectively. ²³⁵U decay chain was divided into ²³⁵U_up, ²³¹Pa, and ²³⁵U_down with ²³⁵U, ²³¹Pa, and ²²⁷Ac as the parent nuclides, respectively. We developed the $\alpha$-ray and $\beta$($\gamma$)-ray background model spectra for each equilibrium section, as shown in Figs. 2 and 3.

⁴⁰K is an environmental isotope with a natural abundance of 0.0117$\%$. The most common disintegration mode is the beta-minus decay to the ground state of ⁴⁰Ca, with a $Q_{\beta^{-}}$ of 1.31 MeV, occurring with a branching ratio of 89.25$\%$. The next dominant mode is the electron capture, which converts ⁴⁰K to the excited state of ⁴⁰Ar, followed by the de-excitation to the ground state with the emission of a 1.46 MeV $\gamma$-ray, occurring with a branching ratio of 10.55$\%$. There are also other modes with branching ratios on the order of less than 0.1$\%$ 40K . Based on the above, the energy spectrum of ⁴⁰K has a photoelectric peak at 1.46 MeV and a continuous beta spectrum up to around 1.31 MeV. Therefore, the radiation emitted from ⁴⁰K can contribute to the background in the 2$\nu$2$\beta$ search because its spectrum is continuous up to around 1.73 MeV. In contrast, ⁴⁰K would not contribute significantly to the background for the search of 0$\nu$2$\beta$, which is expected to produce a peak spectrum at 1.73 MeV.
⁴⁰K is present in both PMT (⁴⁰K_ext.) and GAGG (⁴⁰K_int.). The ⁴⁰K_ext. emits 1.46 MeV $\gamma$-ray, which is detected as background when it enters and deposits the energy in GAGG. Assuming that ⁴⁰K is contained on the photocathode surface of the PMT, we simulated spherically symmetric emissions of 1.46 MeV $\gamma$-rays from a uniform positions on the 500 mm-diameter surface of the light guide model, which is the joint surface with the PMT. For ⁴⁰K_int., we generated stationary ⁴⁰K uniformly in the crystal model and simulated the aforementioned disintegration scheme.
Consequently, we obtained the background model spectra for ⁴⁰K, as shown in Fig. 4. The energy resolution was applied as described in the previous section.

The 2$\beta$ in GAGG can be simulated by generating two electrons from the same position, with their energies governed by theoretical probability distributions. The kinematics of the emitted electrons are described in terms of phase space factors. Reference Phase Space calculated phase space factors for various nuclear candidates across a broad range of mass numbers, from ⁴⁸Ca to ²³⁸U. Reference Phase Space also provides a list of combinations of each electron energy ($\epsilon_{1,2}$) minus rest mass energy ($m_{e}c^{2}$) emitted in the 2$\nu$2$\beta$, denoted as ($\epsilon_{1}-m_{e}c^{2}$, $\epsilon_{2}-m_{e}c^{2}$), along with their corresponding probabilities of occurrence, given at intervals of 0.001 MeV. By summing the probabilities in the list specific to ¹⁶⁰Gd, we produced a discrete cumulative distribution function (CDF) of ($\epsilon_{1}-m_{e}c^{2}$, $\epsilon_{2}-m_{e}c^{2}$). The CDF can then be used to determine ($\epsilon_{1}-m_{e}c^{2}$, $\epsilon_{2}-m_{e}c^{2}$) based on uniform random numbers between 0 and 1, like an inversion method inversion .
The blue histograms in Fig. 5 represent the energy distributions of the summed electrons, $\epsilon_{1}+\epsilon_{2}-2m_{e}c^{2}$ (Left) and each individual electron, $\epsilon_{1}-m_{e}c^{2}$ (Center) and $\epsilon_{2}-m_{e}c^{2}$ (Right), calculated 700,000 times using the CDF. Figure 5 shows that the CDF can determine the energies of the two electrons in accordance with the theoretical distributions represented by the red lines. In GEANT4, we defined two primary electrons emitted with energies $\epsilon_{1}$ and $\epsilon_{2}$ obtained from the CDF, emitted in random directions from a random position in the crystal model. The resulting 2$\nu$2$\beta$ model spectrum is shown on the left in Fig. 6.
The theoretical energy distribution for a single electron of the two electrons emitted in 0$\nu$2$\beta$ is also provided by Ref. Phase Space . In this case, we determined only $\epsilon_{1}$ from the CDF, and calculated $\epsilon_{2}$ by subtracting $\epsilon_{1}$ from the $Q$-value. The 0$\nu$2$\beta$ model spectrum is shown on the right in Fig. 6. It has a peak at the $Q$-value (1.73 MeV) and also has a tail on the low-energy side due to the escape of electrons from the crystal model. The energy resolution was applied as described in Section 2.1.

In this method, we fitted all the pseudo spectra with the background models combined with each 2$\beta$ signal model. We regarded the error in the signal parameter as the upper limits for the number of observed signals ($N_{\rm{obs}}$). We then calculated the lower limits of $T^{2\nu}_{1/2}$ and $T^{0\nu}_{1/2}$ using

In this method, we investigated the $Goodness$-$of$-$Fit$ ($\chi_{2\beta}^{2}$) for the pseudo $\beta$-ray spectrum by increasing signal rate from zero. The amount of U/Th impurities in GAGG was determined from each $\alpha$-ray pseudo spectrum. ⁴⁰K_int. and ⁴⁰K_ext. are fitting parameters associated with each signal rate. The $\chi_{2\beta}^{2}$ is given by

Table 2 summarizes the lower limits of $T^{2\nu}_{1/2}$ and $T^{0\nu}_{1/2}$ in Phase 1, estimated by two different approaches, along with the limits from the previous study. As shown in Table 2, both $T^{2\nu}_{1/2}$ and $T^{0\nu}_{1/2}$ were consistently evaluated for Method 1 and Method 2, and they exceed the limits from the previous study.

The goal of Phase 2 is to either discover 2$\nu$2$\beta$ or constrain the theoretical model for NME, by searching for 2$\beta$ with a sensitivity more than 7.4 $\times$ 10²⁰ years, which is the theoretical prediction Hinohara as described in Section 1. For the 2$\nu$2$\beta$ search, the main radioactive background sources in GAGG are currently 125 $\pm$ 2 mBq/kg of ²³⁸U_up, 14.3 $\pm$ 1.1 mBq/kg of ²³²Th, and 13.0 $\pm$ 1.6 mBq/kg of ⁴⁰K_int. (see Fig. 7(b)). In contrast, the Gd₂O₃, which is the main raw material for GAGG, has been measured to contain $<$ 16.3 mBq/kg, $<$ 0.96 mBq/kg and $<$ 2.7 mBq/kg of ²³⁸U_up, ²³²Th, and ⁴⁰K_int., respectively, as measured by a high purity germanium detector PIKACHU . Other raw materials have also been confirmed to exhibit significantly higher purity than the crystal itself. These discrepancies suggest the contamination from the materials installed in the crystal growing furnace, but they also indicate the potential to develop GAGG with a purity at least an order of magnitude higher than the current level. If such an appropriate purification can be achieved, ⁴⁰K_ext. would become a dominant background in the 2$\nu$2$\beta$ search. We will address this issue through the strategies as described in the next section.
Suppose Phase 2 experiment is conducted for 5 years, using a hundred of GAGG crystals with an order of magnitude lower concentrations of ²³⁸U_up, ²³²Th, and ⁴⁰K_int. compared to those in high-purity GAGG, and utilizing ideal photodetectors without ⁴⁰K_ext.. Under these assumptions, the background spectrum of Phase 2 is expected to be as shown in Fig. 11, which represents one of the $\beta$-ray pseudo spectrum generated in the same way as described in Section 3. By the Method 2 analysis for the 10,000 data-set of the Phase 2 pseudo spectra, we obtained 1.91 (2.84) $\times$ 10²¹ years as the lower limit of $T^{2\nu}_{1/2}$ at 90$\%$ (68$\%$) C.L., providing sufficient sensitivity to meet the objectives of Phase 2.

Crystal purification, including a reassessment of the crystal growth environment and optimization of growth conditions, will remain as a long-term challenge. However, the reduction of ⁴⁰K_ext. can effectively contribute to the indirect improvement of the 2$\nu$2$\beta$ sensitivity, and the strategies for achieving this are relatively straightforward.
⁴⁰K_int. will be a critical background source in the 2$\nu$2$\beta$ search because the beta-minus decay events of ⁴⁰K_int. with a $Q_{\beta^{-}}$ of 1.31 MeV are distributed within the 2$\nu$2$\beta$ energy region, resembling the shape of the 2$\nu$2$\beta$ spectrum. Indeed, their negative correlation can be seen in Fig. 12 (a), which shows the 2$\nu$2$\beta$ rate and the ⁴⁰K_int. radioactivity obtained through the best-fit of all the pseudo data. Therefore, an accurate evaluation of ⁴⁰K_int. is desirable for the precise 2$\nu$2$\beta$ analysis. Additionally, ⁴⁰K_ext. is also negatively correlated with ⁴⁰K_int. due to the similarity of the shape of their model spectra, as shown in Fig. 12 (b), which contributes to the uncertainty of ⁴⁰K_int.. In other words, if the detection rate of ⁴⁰K_ext. can be sufficiently suppressed relative to the radioactivity of ⁴⁰K_int., the ⁴⁰K_int. content can be estimated accurately from the photoelectron peak at 1.46 MeV.

Two main strategies are being considered to reduce ⁴⁰K_ext., the $\gamma$-ray background from ⁴⁰K in the photodetectors. One is to use low radioactivity PMT, while the other is to use multi-pixel photon counter (MPPC) to detect the GAGG scintillation.
Low radioactivity PMTs, such as the R8778 used in the XMASS-I experiment r8778 , have been developed for the dark matter search in low background environments. The R6231-100 was measured to contain 5.58 $\pm$ 0.05 Bq/PMT of ⁴⁰K, as determined by a high purity germanium detector, while the R8778 contains 140 $\pm$ 20 mBq/PMT of ⁴⁰K r8778 , which can suppress the rate of ⁴⁰K_ext. by more than an order of magnitude compared to the current level. However, the quantum efficiency (Q.E.) of R8778 is optimized for about 178 nm wavelength of liquid xenon scintillation Xe . For Q.E. at around 520 nm, corresponding to the wavelength of GAGG scintillation GAGG , the Hamamatsu R6231-04 can be a suitable candidate, as its Bialkali photocathode is sensitive to wavelengths around 300 $\sim$ 650 nm. Furthermore, due to its input window being made of potassium-free glass, the ⁴⁰K content is suppressed to 1.15 $\pm$ 0.07 Bq/PMT, as confirmed by measurement using the high purity germanium detector.
MPPC can also be considered as a photodetector candidate in our experiment, because its peak sensitivity wavelength of it is typically around 450 nm. In fact, GAGG scintillation has been successfully read out using an MPPC array MPPC . The MPPC is composed of fewer materials than the PMT, and in particular, it does not include glass, which typically contains the majority of ⁴⁰K among the PMT components. In this strategy, the MPPC must be arrayed to extend the photosensitive area, thereby covering the cross-sectional area of the GAGG crystal. Additionally, the signal readout circuit must be designed to accommodate the arrayed MPPC.