1 Introduction↩︎

A substantial body of cosmological and astronomical evidence demonstrates that dark matter constitutes a fundamental component of the Universe [1], [2]. The investigation of dark matter is one of the most critical challenges in modern physics. Weakly Interacting Massive Particles (WIMPs) are the most popular dark matter candidates [1]. Numerous experiments have been dedicated to the direct detection of WIMPs, such as XENON [3], PandaX [4], LUX-ZEPLIN [5], DarkSide [6], SuperCDMS [7], EDELWEISS [8], CRESST [9], DAMA/LIBRA [10], and CDEX [11]–[21]. However, WIMPs have not been detected to date. One possible reason is that elastic scattering between WIMPs and nuclei is heavily suppressed. Thus, the Inelastic Dark Matter (iDM) scenario was proposed [22], in which WIMPs scatter with nuclei inelastically.

In inelastic WIMP-nucleus scattering, the WIMP is excited to a higher-energy state (\(\chi^*\)), characterized by a mass splitting \(\delta\) from its ground state (\(\chi\)) [22]. For a given nuclear recoil energy \(E_{\mathrm{nr}}\), there exists a minimal required relative velocity of WIMPs, \[\label{eq:vmin} \begin{align} v_{\min }=\frac{1}{\sqrt{2 E_{\mathrm{nr}} m_{\mathrm{N}}}}\left(\frac{E_{\mathrm{nr}} m_{\mathrm{N}}}{\mu}+\delta\right), \end{align}\tag{1}\] where \(m_{\mathrm{N}}\) is the mass of the target nucleus, and \(\mu\) is the reduced mass of the system. When \(\delta = 0\), the scattering becomes elastic. If \(E_{\mathrm{nr}}\) is too large or too small, \(v_{\min}\) will exceed the maximum velocity of WIMPs that can reach the laboratory, thereby preventing such events from occurring. The maximum and minimum possible values of \(E_{\mathrm{nr}}\) for a certain maximum WIMP velocity \(v_{\max}\) are given by \[\label{eq:Enr} \begin{align} E_{\mathrm{nr,\;max/min}} = \frac{\mu^2 v_{\max}^2}{2m_{\mathrm{N}}} \left( 1 \pm \sqrt{1 - \frac{2\delta}{\mu v_{\max}^2}} \right)^2. \end{align}\tag{2}\] Due to the limits on \(E_{\mathrm{nr}}\), iDM exhibits greater sensitivity to the velocity distribution of dark matter than elastic dark matter.

Moreover, for a given \(\delta\), a minimum WIMP velocity is required for the occurrence of the scattering, \[\label{eq:vgm} \begin{align} v_{\min}^{*} = \sqrt{\frac{2\delta}{\mu}}, \end{align}\tag{3}\] which indicates that if \(\delta\) exceeds a certain threshold, the inelastic scattering will not occur in laboratory experiments due to excessive \(v_{\mathrm{min}}^{*}\).

Analogously to the case of elastic WIMPs, the interaction between inelastic WIMPs and nuclei at initial state (\(A_i\)) could be either spin-independent (SI) [22], [23] or spin-dependent (SD) [24], [25]. This paper focuses exclusively on the SI scenario. In SI scenario, for \(\delta < 1.022\) MeV, the de-excitation of excited-state WIMPs is considered to solely release neutrino-antineutrino pairs [23], which are undetectable in conventional dark matter detectors. Consequently, the nuclear recoil energy is the only observable signature. According to Eq.@eq:eq:vgm , inelastic scattering of Galactic WIMPs with \(\delta \geq 1.022\) MeV is kinematically forbidden in most experimental setups due to the Galactic escape velocity [26] constraint. That is, the physics channel for this analysis is \[\label{eq:chi95A} \begin{gather} \chi + A_i \;\rightarrow\; \chi^{*} + A_f, \\ \chi^{*} \;\rightarrow\; \chi + \nu + \bar{\nu}, \end{gather}\tag{4}\] where the possible nuclear recoil energy \(E_{\mathrm{nr}}\) of the nuclei at their final state (\(A_f\)) are the measureables.

In this study, we place constraints on the inelastic WIMP-nucleon SI interactions with 737.1 kg\(\cdot\)day of data from the CDEX-1B experiment [15] at the China Jinping Underground Laboratory (CJPL) [27]. The devised methodology can be adopted to study a class of iDM models, such as Magnetic Inelastic Dark Matter [24], Effective Field Theory Inelastic Dark Matter [25], Inelastic Dirac Dark Matter [28] and Inelastic Boosted Dark Matter [29], [30].

2 Expected iDM spectra↩︎

The differential nuclear recoil spectrum of inelastic WIMP-nucleus scattering is given by \[\label{eq:dRdEnr} \begin{align} \frac{\mathrm{d}R}{\mathrm{d}E_{\mathrm{nr}}} = \frac{\rho N_{\mathrm{T}}}{m_{\chi}} \int_{v_{\min}}^{v_{\max}} \frac{\mathrm{d}\sigma}{\mathrm{d}E_{\mathrm{nr}}} v f(\vec{v}, t) \mathrm{d}^3 v, \end{align}\tag{5}\] where \(N_{\mathrm{T}}\) is the number of target nuclei per unit effective mass of the detector, \(\rho\) is the local density of WIMPs, \(m_\chi\) is the mass of WIMPs, \(d\sigma/dE\) is the differential cross section of the inelastic scattering, \(f(v)\) is the velocity distribution of WIMPs in the rest frame of the Earth, the lower limit \(v_{\min}\) is formulated by Eq. 1 , the upper limit \(v_{\max}\) is determined by the escape velocity of WIMPs in the Galaxy and the velocity of the Earth [22], [31]. Assuming that the inelastic scattering is spin and energy independent, the differential cross section can be expressed as \[\label{eq:dSigmadE} \begin{align} \frac{d\sigma}{dE_{\mathrm{nr}}} = \frac{m_{\mathrm{N}} \sigma_{\mathrm{n}} }{2 \mu^2 v^2} \cdot (Z \cdot f^{\mathrm{p}} + (A-Z) \cdot f^{\mathrm{n}})^2 F^2_{\mathrm{SI}}(E_{\mathrm{nr}}), \end{align}\tag{6}\] where \(\sigma_{\mathrm{n}}\) is the SI WIMP-nucleon cross section, \(\mu\) is the reduced mass of the WIMP-nucleon system, \(Z\) is the atomic number of the nucleus, \(A\) is the mass number of the nucleus, \(f^{\mathrm{p, n}}\) are the effective WIMP couplings to the proton and neutron, \(F_{\mathrm{SI}}\) is the SI nuclear form factor [22], [31].

High-purity germanium (HPGe) detectors [32], [33], owing to their good energy resolutions and ultra-low energy thresholds, have been applied in dark matter direct detection by CDEX [11]–[21], [34]–[44]. When WIMPs scatter with germanium (Ge) nuclei in the HPGe detector, a portion of the nuclear recoil energy will be converted into detectable ionization energy. This converted energy, called the electron-equivalent energy, is given by \(E_{\mathrm{det}} = Q_{\mathrm{nr}}(E_{\mathrm{nr}}) \cdot E_{\mathrm{nr}}\), where \(Q_{\mathrm{nr}}\) denotes the quenching factor [33]. The differential electron-equivalent spectrum of inelastic WIMP-nucleus scattering in the HPGe detector is given by \[\label{eq:dRdEdet} \begin{align} \frac{dR}{dE_{\mathrm{det}}} = \left( \frac{dQ^{-1}_{\mathrm{nr}}(E_{\mathrm{det}})}{dE_{\mathrm{det}}} E_{\mathrm{det}} + Q^{-1}_{\mathrm{nr}}(E_{\mathrm{det}}) \right) \frac{dR}{dE_{\mathrm{nr}}}, \end{align}\tag{7}\] where \(Q_{\mathrm{nr}}^{-1}\) is the inverse function of the quenching factor function.

We adopt the Standard Halo Model [26]. Accordingly, the WIMP density \(\rho\) is fixed at 0.3 GeV/\((c^2\mathrm{cm}^3)\), the WIMP velocity distribution follows a Maxwellian profile with most probable velocity \(v_0 = 238\) km/s, and the Galactic escape velocity is set to 544 km/s. For the nuclear form factor, the Helm parametrization [45], [46] is employed. The quenching factor \(Q_{\mathrm{nr}}\) is computed using the \(\tt TRIM\) simulation package [47].

Figure 1 shows the expected nuclear recoil and electron-equivalent spectra of inelastic WIMP-nucleus scattering in the HPGe detector, assuming \(m_\chi = 100\) GeV, \(\delta = 100\) keV, and \(\sigma_{\mathrm{n}} = 10^{-40}\;\mathrm{cm}^2\). In the electron-equivalent spectrum, we incorporate the energy resolution of the HPGe detector applied in the CDEX-1B experiment [15].

Figure 2 shows the expected spectra of inelastic WIMP-nucleus scattering in the HPGe detector for different \(\delta\) values, with fixed \(m_\chi\) and \(\sigma_{\mathrm{n}}\). It also shows the expected spectrum for elastic scattering for comparison. From Fig. 2, it can be observed that the distinction between inelastic and elastic scattering on the expected nuclear recoil spectrum primarily manifests as a heavy suppression in the low-energy region of the inelastic scattering spectrum. In addition, Fig. 2 demonstrates that with increasing \(\delta\), the event rate of the entire nuclear recoil spectrum decreases, and the suppressed region expands to higher energies.

Figure 3 shows the expected inelastic scattering spectra in the HPGe detector for different \(m_\chi\) values, with fixed \(\delta\) and \(\sigma_{\mathrm{n}}\). As illustrated in Fig. 3, the nuclear recoil spectrum falls off less sharply as \(m_\chi\) increases.

3 CDEX-1B experiment↩︎

The CDEX experiment utilizes HPGe detectors for direct dark matter detection [11]–[21], [34]–[44] at CJPL [27], located beneath 2400 meters of rock overburden. The CDEX-1B experiment [15] runs one p-type point contact HPGe detector with a 1008 g target mass (fiducial mass of 939 g, after corrections due to a 0.88 \(\pm\) 0.12 mm surface layer) cooled by a cold finger. Additionally, a NaI(Tl) anti-Compton detector is employed to veto background events. With this set up, the experiment demonstrated excellent time stability, low background level, and good energy resolution [15], [18].

The background spectrum of the CDEX-1B experiment, based on the 737.1 kg\(\cdot\)day dataset collected from March 2014 to July 2017 [15], is shown in Fig. 4. The spectrum covers an energy range of 1.5–200 keVee and is processed using anti-Compton veto and bulk-surface event discrimination [15].

4 CDEX-1B background model↩︎

To set constraints on inelastic WIMP dark matter from the CDEX-1B experiment, we construct the background model of CDEX-1B. We first obtain simulated spectra for each background component in the CDEX-1B experiment using \(\tt Geant4\) [48]. Thereafter, we determine the intensities of each component by fitting with maximum likelihood method combined with Markov Chain Monte Carlo (MCMC) algorithm [49]. This process yields an accurate background model of CDEX-1B.

4.1 Background simulation↩︎

The background simulation of CDEX-1B is conducted with \(\tt Geant4\) [48]. In the simulation, we establish a complete CDEX-1B detector model and simulate the decays of radionuclides within detector components.

The structural configuration of the detector comprises several critical components: the Ge crystal, the crystal support structure, the signal pin, the front-end electronics, the vacuum chamber, and the NaI(Tl) anti-Compton detector. The dead layer of the Ge crystal is set to 0.88 mm according to measurement [50].

We account for the following background sources in the simulation: cosmogenic radionuclides in the Ge crystal, radionuclides in surrounding structural materials, and radon progeny located on the outer surface of the vacuum chamber. All simulated radionuclides and their corresponding detector components are listed in Table [tab:bkg95component]. The “Additional \(^{210}\)Pb” in Table [tab:bkg95component] refers to the additional \(^{210}\)Pb contained in the lead materials inside the detector vacuum chamber, independent of the “U Series” and the “Rn Progeny”.

We simulate at least \(10^8\) decay events per background component to ensure statistical validity. The simulation also includes the awgnti-Compton veto of the NaI(Tl) anti-Compton detector [15], bulk-surface event discrimination [50], and the energy resolution of the CDEX-1B detector.

By employing the background simulation, we derive the spectra for the background components of the CDEX-1B experiment.

4.2 Background fitting↩︎

After obtaining the simulated spectra for each background component in CDEX-1B, we fit the experimental spectrum using the simulated spectra to determine their intensities.

We first use the maximum likelihood method implemented with the \(\tt scipy.optimize\) module [51] of Python to fit the spectrum. However, the fitting results obtained through this method remain insufficiently accurate.

Thus, we further employ the \(\tt emcee\) implementation [52] of the MCMC algorithm [49] to conduct higher-accuracy iterative refinements based on the initial optimization results.

To eliminate the influence of discrepancies in peak morphology between the experimental and simulated spectra on the fitting procedure, we merge each peak region defined by the identified peaks in the experimental spectrum into a single composite bins during spectral fitting.

Using the MCMC algorithm [49], we refine the fit to achieve higher accuracy. This enables us to establish an accurate background model for the CDEX-1B experiment. As shown in Fig. 5, the background model shows good agreement with the experimental background spectrum, with the majority of bins exhibiting deviations within the 3 \(\sigma\) error range.

5 Constraints on iDM↩︎

Based on the expected spectra of iDM in the CDEX-1B detector and the background model of the CDEX-1B experiment, we calculate constraints on SI inelastic WIMPs from the CDEX-1B experiment.

5.1 Calculation method↩︎

Following the methodology established during the construction of the background model, we employ a similar method based on the maximum likelihood fitting method and the MCMC algorithm [49] to calculate constraints.

Firstly, we fit the experimental background spectrum of CDEX-1B using simulated spectra of each background component combined with the expected spectrum of inelastic WIMPs with given mass and splitting energy. The fitting utilizes the maximum likelihood method and the \(\tt scipy.optimize\) module [51]. The fit parameters include the intensities of background components and the WIMP signal, with initial values for the background intensities taken from the background model.

Subsequently, we apply the \(\tt emcee\) implementation [52] of the MCMC algorithm to perform Bayesian parameter estimation, sampling the posterior distribution starting from the initial optimization results.

In both fitting stages, we set the lower limit of the WIMP signal intensity to zero, meaning that negative intensities are excluded. Additionally, following the same methodology employed during the construction of the background model, we merge each peak region defined by the identified peaks in the experimental spectrum into a single composite bins during spectral fitting.

By adopting the MCMC algorithm, we obtain the best-fit result and the posterior distribution of the WIMP signal intensity. Based on the posterior distribution of the WIMP signal intensity, we derive the posterior distribution of the corresponding \(\sigma_{\mathrm{n}}\). The 90% quantile of the posterior distribution of \(\sigma_{\mathrm{n}}\) represents the 90% confidence level (C.L.) upper limit (one-sided) on it.

Figure 6 illustrates the background and inelastic WIMP spectra corresponding to the best-fit result, the inelastic WIMP spectrum corresponding to the 90% C.L. upper limit, and the posterior distribution of \(\sigma_{\mathrm{n}}\) at \(m_\chi = 100\) GeV, \(\delta = 100\) keV.

5.2 Results↩︎

Using the above method, we calculate constraints on inelastic WIMPs for various \(m_\chi\) and \(\delta\) values. The low 1.5 keVee analysis threshold of CDEX-1B allows the studies of the expanded parameter space of \(\delta \ge 2\) keV and \(m_\chi \ge 10\) GeV. Moreover, as can be seen from Eq. 3 , constrained by the Galactic escape velocity of WIMPs [26], there exists an upper limit to the detectable \(\delta\) for any given \(m_\chi\) values.

Figure 7 shows the 90% C.L. upper limits on \(\sigma_{\mathrm{n}}\) from the CDEX-1B experiment for \(m_\chi = 10\), 25, 50, 100, 250, and 500 GeV. For \(m_\chi = 100\), 250, and 500 GeV, we also present the 90% C.L. upper limits from CDMS[53], [54], CRESST [54], [55], and ZEPLIN-I [54], [56], along with the allowed regions from DAMA/LIBRA [54], [57], [58] for comparison. At \(m_\chi = 250\) GeV, the CDEX-1B upper limits fully exclude the \(\chi^2 < 4\) allowed regions from DAMA/LIBRA for \(\delta < 30\) keV. At \(m_\chi = 500\) GeV, the CDEX-1B upper limits fully exclude the \(\chi^2 < 9\) allowed region from DAMA/LIBRA for \(\delta < 50\) keV.

6 Conclusions and discussions↩︎

In this work, we present the SI inelastic WIMP dark matter search results from the CDEX-1B experiment [15]. By establishing an accurate background model of the CDEX-1B experiment, we calculate constraints on inelastic WIMPs using the CDEX-1B dataset with an exposure of 737.1 kg\(\cdot\)day [15]. The \(90\%\) C.L. upper limits on the SI WIMP-nucleon cross section \(\sigma_{\mathrm{n}}\) from CDEX-1B fully exclude the \(\chi^2 < 4\) allowed regions for \(\delta <\) 30 keV at \(m_\chi =\) 250 GeV and the \(\chi^2 < 9\) allowed region for \(\delta <\) 50 keV at \(m_\chi\) = 500 GeV from DAMA/LIBRA [54], [57], [58].

A method is developed in this paper to establish the CDEX-1B background model and to place constraints on iDM. The analysis procedures can be adopted to probe various variants of iDM models [24], [25], [28]–[30].

Future studies could explore the detection of decay products from excited iDM [30], [59]. Additionally, since iDM is more sensitive to the velocity distribution of dark matter compared to elastic dark matter, better results may be achieved by using the modulation method for analysis [18], [22], [60].

The next-generation of the CDEX experiment, CDEX-50 [20], is being prepared. In the CDEX-50 experiment, an upgraded detector array consisting of 50 HPGe detectors with a target mass of 50 kg will be deployed. The CDEX-50 experiment has a target exposure of 150 kg\(\cdot\)year, and its background level will be reduced to approximately \(2 \times 10^{-4}\) \(\mathrm{counts \cdot kg^{-1} \cdot keVee^{-1} \cdot day^{-1}}\) (cpkkd) at 20 keVee [20], which is \(5 \times 10^3\) times lower than that of CDEX-1B. Therefore, the results of the CDEX-50 experiment are expected to improve limits on \(\sigma_{\mathrm{n}}\) by 4 orders of magnitude compared to our current results.

This work was supported by the National Key Research and Development Program of China (Grants No. 2017YFA0402200 and No. 2022YFA1605000) and the National Natural Science Foundation of China (Grants No. 12322511, No. 12175112, No. 12005111, and No. 11725522). We acknowledge the Center of High Performance Computing, Tsinghua University, for providing the facility support. We would like to thank CJPL and its staff for hosting and supporting the CDEX project. CJPL is jointly operated by Tsinghua University and Yalong River Hydropower Development Company.

References↩︎