Nonproportionality of NaI(Tl) Scintillation Detector for Dark Matter Search Experiments

We present a comprehensive study of the non-proportionality of NaI(Tl) scintillation detectors within the context of dark matter search experiments. Our investigation, which integrates COSINE-100 data with supplementary γ spectroscopy, measures light yields across diverse energy levels from full-energy γ peaks produced by the decays of various isotopes. These γ peaks of interest were produced by decays supported by both long and short-lived isotopes. Analyzing peaks from decays supported only by short-lived isotopes presented a unique challenge due to their limited statistics and overlapping energies, which was overcome by long-term data collection and a time-dependent analysis. A key achievement is the direct measurement of the 0.87 keV light yield, resulting from the cascade following electron

In recent years, the dark matter (DM) community has shown increasing interest in the nPR of NaI(Tl), particularly in the low-energy region.It is driven by experiments searching for DM using NaI(Tl) detectors, inspired by the observations made by DAMA/LIBRA [32][33][34][35].DAMA/LIBRA claimed to have detected an annual modulation signal compatible with DM interaction below 6 keV, which has motivated other experiments to analyze low-energy scintillation events in an attempt to replicate their results [36,37].Moreover, the search for low-mass DM signal in the low-energy region has intensified [38,39], with consideration of interactions such as weakly interacting massive particle (WIMP) scattering on the nucleus [40][41][42][43][44][45], the Migdal effect [46][47][48], and DMelectron scattering [49][50][51][52][53][54][55].Another phenomenon that relies on the analysis of low-energy scintillation is coherent elastic neutrino-nucleus scattering (CEνNS).Ongoing efforts aim to observe CEνNS of reactor neutrinos using NaI(Tl) detectors [56].Additionally, a feasibility test of future NaI(Tl) detectors for CEνNS observations of solar and supernova neutrinos has been conducted [57].These endeavors crucially require understanding extremely low-energy events.
While substantial efforts have been dedicated to lowering the energy threshold through hardware and software upgrades [36,58,59], there has been a relative neglect of the study on nPR of scintillation response.However, misconceptions regarding scintillation response could lead to incorrect interpretations of DM and ν signals.It is also evident in the recent measurements and interpretations of nuclear recoil quenching factors, where taking nPR into account results in an energy scale difference of about 25% at 1 keV electron equivalent [60].This is particularly relevant in the low-energy region, where light yield rapidly changes and can introduce significant biases.Unfortunately, previously measured nPR curves have not adequately covered the low-energy region, resulting in a scarcity of reported light-yield measurements at these energies, with only a few indirect or obsolete values available [3,10].
To address the aforementioned challenges, this paper presents the results of our γ-ray calibration of NaI(Tl) crystals used in the COSINE-100 experiment.We utilized various mono-energetic peaks from internally contaminated isotopes, covering an energy range from 0.87 to 88 keV.Given the long-term nature of the DM search experiment, we extracted these energy peaks by simultaneously examining their decay times and peak positions.To complement the data, we also conducted γ spectroscopy using another NaI(Tl) crystal from the same manufacturer.
In Sec. 2, we describe the experimental setup, encompassing the COSINE-100 experiment, as well as the additional γ spectroscopy.COSINE-100 data, which forms the basis of this study, is summarized in Sec. 3, highlighting the observations and fitting procedures.Sec. 4 focuses on the γ spectroscopy methodology, outlining the measurements and analysis techniques to supplement COSINE-100 data.In Sec. 5, we consolidate the results and present the measured nPR.Additionally, we provide the resolution and validate it using a waveform simulation.Finally, in Sec.6, we conclude by summarizing the key findings of our research and emphasizing the implication for future research.

Experimental Setup
The COSINE-100 experiment employed eight high-purity NaI(Tl) crystalsdeveloped in collaboration with Alpha Spectra Inc.These crystals collectively served as the target for DM direct detection.They had a total mass of 106 kg and were equipped with one 3-inch photomultiplier tube (PMT) at each of the two ends, enabling the detection of scintillation light.Each crystal was housed within a copper encapsulation.The setup featured eight such encapsulations submerged in a 2,200 L volume of LAB-based liquid scintillator (LS), which served as an active veto system [61].An acrylic box contained the LS, with eighteen PMTs attached to its walls for scintillation light detection, respectively.Events were categorized as multiple-hit or single-hit events based on whether particles left signals in other detectors (LS or crystal).The acrylic box is surrounded by copper and lead, which provide passive shielding layers, while an outermost layer consists of plastic scintillator panels used as a muon veto [62].
Signals collected by the crystal PMTs were recorded, as waveforms, through a data acquisition (DAQ) system.The DAQ system comprised pre-amplifiers, analog-to-digital converters (ADC), and a trigger control board (TCB), with detailed specifications outlined in Sec.7 of Ref. [63].The charge deposition of an event, serving as an indicator of energy, is determined through the integration of the waveform over a 5 µs time window following noise removal.Data taking for physics purposes commenced on October 20, 2016, at the Yangyang Underground Laboratory [63], and this study utilized data collected through June 9, 2022.Notably, three crystals were excluded from the analysis due to their low light yields or high noise levels [42,44].
Despite the high purity of these NaI(Tl) crystals, radioactive decay signals were observed inside them [64].Specifically, in the region of interest of this study, three long-lived isotopes -210 Pb, 22 Na, and 40 K -were clearly identified as contamination inside the crystals.Moreover, several short-lived isotopes 1 were generated by cosmic rays prior to installation, including 109 Cd, 113 Sn, 125 I, 121m Te, and 127m Te [65].Signal rates from these isotopes were characterized by the isotopes' decay times and effectively modeled using Monte Carlo simulations.This paper reanalyzes the energy peaks associated with these isotopes, which serve as calibration points.
Apart from the crystals used in COSINE-100, a smallersized NaI(Tl) crystal from the same manufacturer was utilized for a γ spectroscopy experiment.The crystal had a cylindrical shape with a diameter of 3 inches and a height of 7 inches.It was enveloped in 2-3 mm of PTFE, enclosed by 2 mm of aluminum, and fitted with a quartz window.It featured a 3-inch PMT identical to those used in COSINE-100.The crystal was further shielded using 5-cm-thick lead bricks and installed at the IBS HQ ground laboratory in Daejeon, Korea.The DAQ system for this crystal also included a pre-amplifier, an ADC module, and a TCB, all with the same specifications as those in COSINE-100.Notably, a more stringent trigger condition was applied to reduce the high trigger rate arising from noise and background events at the lowenergy region.This setup facilitated measurements of energy responses to external γ sources, such as 241 Am, 133 Ba, 109 Cd, and 137 Cs.Importantly, the peaks used for analysis remained independent of the trigger condition.

Extracting Internal Peaks
From COSINE-100 data, there are peaks resulting from the decay of isotopes contaminated in the crystals, which were used for calibration.Throughout the entire experimental period, three peaks from long-lived isotopes -49 keV from210 Pb, 0.87 keV from 22 Na, and 3.2 keV from 40 Kremained consistently visible.These peaks were extracted by modeling their shapes and the background components around them from the charge distributions stacked over the entire experiment duration.Short-lived isotopes also produced additional charge peaks that were visible only in the early data.To extract them, a time-dependent model for charge distributions was developed in accordance with the decay of each isotope.The following sections outline the formulation and results of the modeling.

Long-Lived Isotopes
The charge distribution of single-hit events in COSINE-100 revealed a distinct peak from 210 Pb.It decays to 210 Bi through β decay, where 84% of them are excited at 46.5 keV.The subsequent deexcitation of this daughter nucleus emits 1 Precisely to say, short-lived isotopes without long-lived supporting parent isotopes.an energy of 46.5 keV, while the energy of β ray is widely distributed, having an average of 4.2 keV.Due to the asymmetric shape of the β spectrum and the resolution-smearing effect, the peak position becomes approximately 49 keV.
To effectively represent this composition, the peak is modeled with the function f Pb (q), defined as f Pb (q) =Q (q; A, q min , q 0 ) * N q; µ = 0, σ 2 where q is the integrated charge from the signal waveform, while all the other parameters are free fitting parameters.In this equation, the first part (Q * N) represents the 49 keV peak from 210 Pb, where Q is the quadratic distribution, defined as where A is the normalization factor, q 0 is the vertex position, and q min is the minimum range of the distribution.In the context of the decay process of 210 Pb, q min represents the charge deposition associated with the energy of the γ ray, while q 0 corresponds to the charge deposition at the maximum total energy.The symbol * denotes the convolution operator, and N is the normal distribution centered at 0 with its standard deviation given by the parameter σ .The convolution of this normal distribution mimics the resolution-smearing effect in the charge distribution.
For the background distribution around this peak, we adopted a sum of a linear component (Cq + D) and an asym-metric Gaussian distribution AG defined as where B is the normalization constant, µ is the peak position, and σ l and σ r are the standard deviations impacting the left and right tails.The distributions accurately represented the background components, as confirmed by the simulated distributions created using the GEANT4 Monte Carlo simulation toolkit [66].The fitting method produced a bias of less than 0.1% in the peak position, which suggests that it is reliable enough to handle other subdominant radioactive isotopes like 121m Te or 129 I.The data distribution and the fitted model are depicted in Fig. 1.
The position of the 210 Pb peak was determined as the position where the peak component (Q * N) is maximized.To interpret the position as the light amount for the energy, we assign a 1% systematic uncertainty to account for the unknown spatial distribution of 210 Pb within the crystal and for the poorly known spatial dependence of the detector response.The total uncertainty is visually depicted as a yellow band around the vertical line representing the peak position in Fig. 1.The stability of fitting was improved by assigning the peak position based on its maximum value.However, this also caused minor fluctuations in the true peak energies among the detectors due to their resolution difference.The bias term was determined to be about 0.5% and was considered when calculating the light yield.
In the charge distribution of multiple-hit events, we observed two distinct peaks associated with isotopes from internal contamination, specifically 22 Na and 40 K. Approximately 10% of 22 Na decays to the 1275 keV level of 22 Ne through electron capture (EC), while a similar proportion of 40 K decays via EC to the 1461 keV level of 40 Ar [67,68].The daughter atoms undergo cascade processes, leaving characteristic peaks at their respective K-shell edges within the crystal where 22 Na and 40 K originally resided.Concurrently, high-energy γ rays are emitted due to the deexcitation of these daughter nuclei.These energetic γ rays can escape the crystal, potentially triggering the neighboring crystals or the LS veto system, leading to coincidences.Consequently, we clearly identified two low-energy peaks in the charge distribution of multiple-hit events, corresponding to 0.87 keV and 3.2 keV.These values represent the K-shell binding energies of 22 Ne and 40 Ar [69].The charge distribution is shown in Fig. 2, where the observed peak positions and shapes align with our expectations.Fig. 2 The charge distribution of multiple hit events at the low-energy region, as recorded by a NaI(Tl) crystal in COSINE-100 (indicated by the black crosses).The solid orange curve corresponds to the fitted model.Additionally, the dashed yellow curve represents the scaled Poisson distribution for the 0.87 keV from 22 Na, while the dashed pink curve illustrates the normal distribution for the 3.2 keV from 40 K.The vertical lines and bands depict their peak positions and associated errors.The dotted blue curve represents the linear-shape model for the background components.
For the 0.87 keV peak, a scaled Poisson distribution was used for modeling.The distribution is represented as where A is the normalization constant, λ is the expected number of photo-electrons (N pe ), and s denotes the scaling factor originating from the PMTs.The 3.2 keV peak produced sufficient photo-electrons such that we modeled it using a normal distribution.The background components were attributed to Compton scattering of high-energy external γ rays, resulting in relatively smooth shapes.To account for the uncertainty in the background distribution's shape, flat or linear-shaped distributions were adopted as background models, and the difference between results from these models was considered as a systematic uncertainty.All the parameters denoting the normalization, position, and width of the peaks, and the shape of background distribution are free in the fitting process.The charge distribution, fitted model, and the extracted peak positions are shown in Fig. 2, and the peak positions are summarized in Table 1.

Short-Lived Isotopes
Short-lived isotopes, such as 109 Cd, 113 Sn, 125 I, 121m Te, and 127m Te, were initially activated by cosmic rays while the crystals were exposed on the ground and disappeared quickly due to their fast decay times.Their decay processes emit γ rays and X-rays at various energy levels, making them valuable for this study.Unfortunately, their energy levels overlap, and their activities are scarce, posing challenges for extracting their peaks using the method outlined in Sec.3.1.
To overcome this issue, we leveraged constraints based on their decay characteristics.
To observe the decays of these isotopes, we accumulated charge distributions over 475 days from the beginning of data acquisition, dividing it into 19 periods, each spanning 25 days.Fig. 3 shows the distributions for three periods: 1 st , 5 th , and 9 th .While decaying and constant components are visually identifiable, determining the position of each peak directly from the spectral shape is challenging.We built a model to extract the peaks from these spectra, with consideration of the decay characteristics of the short-lived isotopes.
The time-dependent charge distributions for each shortlived isotope i can be expressed as where t denotes time, A i and τ i are the initial activity and lifetime of isotope i, and f i j is the branching fraction for the j-th peak of isotope i.Each peak is modeled as a normal distribution N(q) with peak position µ i j and resolution σ i j .
For the 210 Pb component, the normal distribution in Eq. 5 was replaced with Q * N, as used in Sec.3.1.
The charge distribution accumulated between times t s and t e can be modeled as where s(t) denotes the detector operational state (1 or 0) at time t.The term C(q) represents the time-independent component in the spectrum, consisting of isotopes with decay times orders of magnitude longer than the experiment's duration.Instead of constructing a complex analytic model for the uninteresting C(q) term, we attempted to eliminate its explicit dependence.Since the term is common to all time periods, subtracting a charge distribution from another period with proper time-normalization could effectively remove the C(q) dependency.Fig. 3 The charge distributions from the initial COSINE-100 data and the extracted peaks attributed to short-lived isotopes.The data points, indicated by black crosses, represent measurements from each data collection period and have been simultaneously fitted by using Eq. 8.The time-normalized late distribution L(q) is depicted as gray histograms in each data collection period.Each isotope component is visualized as a stacked histogram in its own color.
For the subtraction, a late distribution L(q) was accumulated, expressed as L (q) = S (q, t 0 , t 1 ) for t 0 and t 1 defined as 650 days and 1,150 days after the start of data acquisition.Due to its long time period, L(q) has a sufficiently small error to be treated in the model as a known distribution rather than fitted data.Replacing the C(q) term with L(q) from Eq. 6 yields S (q, t s , t e ; {A i , ⃗ µ i , ⃗ σ i }) which depends only on the set of parameters related to the peaks.Consequently, the number of events for the k-th 25days-period spectrum in the l-th charge bin is expected as S kl = S (q l , (k − 1)∆ T, k∆ T ), where q l is the central value of the l-th charge bin, and ∆ T is the period interval, 25 days.
To extract peaks from the short-lived isotopes, we conducted a simultaneous fit on 19 distributions using a binned likelihood, where O kl is the observed number of events in the k-th period spectrum at the l-th charge bin.The initial activity A i , peak positions µ i j , and resolutions σ i j were treated as free parameters, while those for 210 Pb were kept constant at the values obtained in Sec.3.1.Constraints for the initial activities, µ A i and σ A i were obtained from our previous measurement [65].As seen in Fig. 3, the extracted peaks effectively account for the charge distributions in different time domains simultaneously.
For verification, we applied the same fitting procedure to simulated distributions, allowing us to obtain the true mean energy value for each peak.Through this process, light yields were obtained for several points ranging from 25 to 88 keV, as summarized in Table 1.Crystal 2, cooled underground for about 2.75 years before the experiment started, had the lowest amount of short-lived cosmogenically induced isotopes [65].It was the longest cooling time among the crystals, resulting in larger errors and missing elements in Table 1 for Crystal 2. The true energy estimation from the simulated distribution for Crystal 2 was also unstable due to the same reason, and these uncertainties were considered in the nPR analysis.Moreover, the nPR analysis excluded the peaks from isotopes 121m Te and 127m Te because their amounts were small and their peak positions overlapped with stronger isotopes.

Gamma Spectroscopy
An independent experiment was conducted for γ spectroscopy to acquire complementary data and broaden the energy range of our analysis.Charge distributions were measured using four γ sources, 241 Am, 133 Ba, 109 Cd, and 137 Cs.These sources were placed near the crystal and aligned at its center during data taking.The energies of the γ rays and X-rays emitted from these sources, along with their exposure times, are listed in Table 2. Data without any sources was also collected over an extended period to obtain the background distributions before and after acquiring data with the sources.This allowed us to verify the stability of the PMT gain and background distribution with time.The pure peak for each source was identified by subtracting the background distribution, as illustrated in Fig 4.
The peak positions were determined by fitting the charge distributions.The charge distribution for isotope i can be modeled as where N i j represents the normalization constant of the j-th peak, and CB(q) denotes a Crystal Ball function commonly used to model peaks measured by scintillation crystals with energy loss processes [70].The function includes parameters for peak position µ i j and resolution σ i j , which are free during the fitting process.The α i j and n i j parameters account for the threshold and the exponent of the power-law low-end tail part in the distribution, respectively.Flat and linear background components were added to the 133 Ba and 137 Cs distributions, respectively, to model Compton scattering of high-energy γ rays.
The model was fitted to data by minimizing the following chi-squared: In this equation, D il and B l are the measured number of events in i-th source and background data, respectively, at the l-th charge bin, centered at q l .Exposure times for the  i-th source and background data are represented by t i and t b , respectively.The extracted peak positions are displayed in Fig. 4 as vertical lines with their uncertainties.The same fitting procedure was conducted with simulated distributions, as done in Sec. 3.All dimensions of the detector geometry were implemented based on measurements, except for the thickness of the PTFE wrapped inside the aluminum case.To account for the uncertainty, two different conditions were simulated with this thicknesses set to 2 mm and 3 mm.The simulation results were fitted separately, and the difference in peak positions was assigned to the systematic error of the true energy.One can see those errors and the extracted peak positions in Table 2.

Scintillation Response of NaI(Tl) Crystal
This section presents the nPR of NaI(Tl) crystals based on the measurements described in the previous sections.We also compare our results with previous references.In addition, the energy resolution for each peak was measured and vali-dated using a waveform simulation developed, considering the characteristics of the NaI(Tl) scintillator [71].

Photon Response Measurements
The relative light yield obtained via peak extraction is depicted in Fig. 5.In the previous sections, we measured numerous charges of peaks from the data and estimated the true energy values from simulated spectra.The charge was divided by the true energy for each peak.To emphasize nonproportional behavior over the absolute values, the ratio of charge to energy was normalized by a constant for each crystal, termed the relative light yield.Normalization was chosen to set the value to unity at 49 keV for all crystals.The error bars include systematic errors resulting from the radioactive source's position and the peak extraction method.These measurements exhibited consistency across the crystals and the measuring techniques within the associated uncertainties.The well-known K-shell dip of NaI(Tl) at 33 keV is also identifiable.
The data was fitted using an empirical function developed based on the Payne's model [15][16][17][18][19], as shown by the orange curve in Fig. 5.An expectation arises for the presence of the L-shell dip structure around 4.5 keV, although there is an insufficient number of data points below 20 keV to constrain the detailed shape of the dip.We have checked the consistency of our curve with measurements from two references comparing quantities, such as the light-yield ratio between 4.5 and 20 keV, the ratio between 0.87 and 20 keV, and the slope around 10 keV.These were found to be consistent with previous reports [3,10], as summarized in Table 3.Despite this consistency check, we acknowledge that our curve between 4.5 and 20 keV is less convincing, resulting in a relatively wider error band, with its width dependent on the selected model parameterization.The measurement of light yield in this energy region is left for future investigation.
To evaluate the nPR in the high-energy region above 200 keV, we examined its consistency with the data reported

Relative Light Yield
Reference Data Fig. 6 The relative light yield (normalized at 662 keV) for the highenergy γ peaks is represented by black dots with error bars.The orange dots and solid orange curve represent the nPR measured from Ref. [7], where our data agrees well with the reference.
in Ref. [7].We measured the light yields by fitting the prominent γ peaks from the high-energy spectrum, including 295 keV from 214 Pb, 511 and 1274 keV from 22 Na, 609 and 1764 keV from 214 Bi, 1461 keV from 40 K, and 2615 keV from 208 Tl.
The measured light yields demonstrate agreement with the nPR curve reported in Ref. [7], as illustrated in Fig. 6.

Resolution and Waveform Simulation
Once a peak is extracted, the energy resolution σ i j is naturally obtained along with the peak position µ i j .As shown in the dashed gray curve in Fig. 7, the characteristic 1/ √ E-like shape, originating from the Poisson fluctuations in N pe for the PMT, is clearly visible.Also resolution effects due to the detector geometry and the distribution of the number of amplified electrons in the PMT are observed.Best Fit χ 2 /dof = 5.36 / (9-2) √ χ / q Compo e t Fig. 7 The relative resolution as a function of incident γ-ray energy for a NaI(Tl) crystal.The black dots with error bars represent the measurements, and the solid orange curve represents the model of the resolution curve, as described in Eq. 11, with an error band at 68% confidence level.The solid pink curve illustrates the a/q component from the resolution model, which closely matches the waveform simulation (solid blue curve).The dashed gray curve represents the Poisson fluctuation in N pe .
To account for these effects, the relative resolution as a function of charge was modeled using the following equation, Here, a and b are fitting parameters, where a scales the Poisson resolution contribution and b typically arises from nonuniformity of detector response, which is not included in the simulation.In Fig. 7, the orange line represents the model fitted to the measurements, with the uncertainty shown as a band.The resolution σ (q) is transformed to σ (E) via the nPR curve shown in Fig. 5, leading to the emergence of non-analytic points at the 4 keV L-shell edge and the 33 keV K-shell edge.The pink line represents the function when b = 0. Given the focus on DM detection experiments, where the region of interest typically lies below 6 keV, it is crucial to investigate the resolution in the sub-keV range, particularly for future upgrades.Since the effect of electron amplification in PMTs and DAQ system is most dominant at low energies, except for that due to the Poisson fluctuations with N pe , a waveform simulation developed for the COSINE-100 experiment was utilized to validate the resolution in this region [71].This simulation incorporates factors such as rise and decay times of NaI(Tl) scintillation, fluctuations in the single photoelectron shape, DAQ electronics, and charge accumulation processes.Therefore, by inputting N pe that takes into account Poisson fluctuations, we can estimate the resolution with these effects included.
The blue line in Fig. 7 shows that the simulation successfully replicates the data.Despite not accounting for the intricate details of detector geometry and other mechanisms, the simulation exhibits good agreement with the measurements at the low-energy region.In addition, the agreement between the simulation and the pink curve, which is the measurement with the constant-term factor excluded, confirms the possibility of extrapolating the resolution curve to the low-energy direction.These findings hold particular significance for scintillator experiments focusing on low-energy ranges.

Summary and Conclusion
In conclusion, this study has thoroughly examined the nPR characteristics of NaI(Tl) scintillators within the context of DM search experiments.Through a comprehensive calibration process using COSINE-100 data, complemented by γ spectroscopy, we have gained valuable insights into the behavior of NaI(Tl) crystals.
The measurement of light yields at various energies was successfully achieved by analyzing internal peaks originating from both long-lived and short-lived isotopes.Extraction of peaks from long-lived isotopes was easily accomplished due to stable experimental data acquisition.A significant milestone was reached by directly measuring the light yield at 0.87 keV, generated by the cascade following 22 Na EC decay, facilitated by the active veto system.This direct measurement serves as validation of an indirect measurement obtained by analyzing X-rays near the Iodine K-shell using synchrotron radiation [10].Overcoming the challenge of extracting multiple overlaid peaks from short-lived isotopes was made possible through long-term data collection and time-dependent analysis.
In the nPR curve measured here, the distinctive K-shell dip structure was observed.The nPR curve was fitted using an empirical curve.The resulting fitted nPR curve exhibits consistency with previously reported measurements, particularly in the energy range below 20 keV.However, the lack of data points around the 10 keV region limits the precision of the calibration curve, leaving the error band relatively wider.This aspect remains a focus for future investigation.
The resolution analysis of NaI(Tl) scintillators revealed a characteristic resolution curve that follows the 1/ N pe trend, as expected for Poisson fluctuations in N pe .Waveform simulation provided compelling confirmation of the resolution at low energies, allowing for the extension of the resolution curve to lower energy ranges.This enhancement in our ability to detect low-energy DM recoil signals is very important.
The findings of this study hold significant implications for background modeling in NaI(Tl) experiments.Precise calibration and characterization of the nPR effect are essential for accurately modeling the background spectra resulting from radioactive interactions involving electrons and photons.To facilitate more detailed simulations, the integration of the nonproportional light production of NaI(Tl) into the GEANT4 toolkit is under development, building on similar efforts reported previously [72,73].
Furthermore, it is imperative to consider nPR calibration when interpreting signals that may come from DM interactions with electrons or photons.As the scientific community explores diverse scenarios involving boosted DM and light DM, the significance of these interactions is poised to grow.This study lays the foundation for more robust and accurate DM search experiments in the future.

Fig. 1
Fig.1The charge distribution around the 49 keV peak of 210 Pb, as measured by a NaI(Tl) crystal in COSINE-100 (depicted by the black crosses).The solid orange curve represents the fitted model described in Eq. 1.The 210 Pb peak is drawn as the dashed yellow curve.The vertical yellow line indicates the position of the 210 Pb peak, along with its associated error.The dotted blue curve represents the summed spectrum of background components.

Fig. 4
Fig. 4 Charge distributions of external sources after subtraction of the background distribution, recorded by the sample crystal for γ spectroscopy.Black crosses represent the data, and orange curves depict the fitting results used to extract the peak positions.The peak positions extracted are shown as vertical lines with error bands.

Fig. 5
Fig.5The relative light yield as a function of incident γ-ray energy, normalized to the value at 49 keV.The measurements are represented by colored dots with error bars, and the fitted function is drawn as the solid orange curve.The translucent orange band represents the uncertainty from the modeling at 68% confidence level.

Table 1
List of peak positions extracted from COSINE-100 data.

Table 2
List of sources used for γ spectroscopy, with measured peak positions.The energy was obtained by modeling the simulated spectra.

Table 3
Comparison of γ nPR measurements.LY γ (E) is the light yield for a γ ray with incident energy E, and the relative light yield LY γ, 49 (E) ≡ LY γ (E) /LY γ (49 keV).