Strong constraints from COSINE-100 on the DAMA dark matter results using the same sodium iodide target.

We present new constraints on dark matter interactions using 1.7 years of COSINE-100 data. The COSINE-100 experiment, consisting of 106 kg of tallium-doped sodium iodide [NaI(Tl)] target material, is aimed to test DAMA’s claim of dark matter observation using the same NaI(Tl) detectors. Improved event selection requirements, a more precise understanding of the detector background, and the use of a larger dataset considerably enhance the COSINE-100 sensitivity for dark matter detection. No signal consistent with the dark matter interaction is identified and rules out model-dependent dark matter interpretations of the DAMA signals in the specific context of standard halo model with the same NaI(Tl) target for various interaction hypotheses.

INTRODUCTION

Astronomical observations continue to indicate that the Universe is made mostly of nonluminous and invisible dark matter (, ). Several types of new fundamental particles have been proposed as candidates for the dark matter () such as weakly interacting massive particles (WIMPs) (, ), but no definitive signal has been seen despite concerted efforts by many collaborations (). One exception is the much-debated claim by the DAMA collaboration of a statistically significant annual modulation in the event rate of their experiment (–) with a period and phase consistent with that expected from WIMP dark matter (–). This is controversial because if it is interpreted as a signature for WIMP interactions, it conflicts with other direct search experiments (–) that report null signals in the regions of parameter space that are allowed by the DAMA signals.

Several groups have been working to develop new experiments with the aim of reproducing or refuting DAMA’s results using the same NaI(Tl) medium (–). The COSINE-100 experiment is one of them and that is currently operating with 106 kg of low-background sodium iodide crystals at the Yangyang Underground Laboratory (, ). An analysis of the initial 59.5 days of COSINE-100 data showed that the annual modulation signal reported by DAMA is inconsistent with the explanation using SI interaction between WIMPs and sodium or iodine nuclei in the context of the standard halo model (SHM) (). However, this first result left open interpretations using certain alternative dark matter models (, ), dark matter halo distributions (), and detector responses (, ) that could allow room for consistency between DAMA and COSINE-100. Model-independent searches of an annual modulation signal using 1.7-years data were also reported but were still not sensitive enough to conclusively challenge the DAMA observation (). Here, we present results from an analysis of 1.7 years of COSINE-100 data with improved event selection requirements and an energy threshold that has been reduced from 2 to 1 keVee, where keVee is electron-equivalent energy (). We find an order-of-magnitude improvement in sensitivity, sufficient to strongly constrain these alternative scenarios, as well as to further strengthen the previously observed inconsistency with the WIMP-nucleon spin-independent interaction hypothesis ().

RESULTS

Experiment

COSINE-100 is located at the Yangyang Underground Laboratory in South Korea with a 700-m rock overburden (, ). The experiment consists of eight low-background thallium-doped sodium iodide [NaI(Tl)] crystals arranged in a 4 × 2 array with a total target mass of 106 kg. The array is immersed in 2200 liters of liquid scintillator used to identify events induced by radioactive backgrounds inside or outside the crystals (). The liquid scintillator is surrounded by copper, lead, and plastic scintillators to reduce the background contribution from external radiation as well as to tag and reject events associated with cosmic-ray muons (). Each NaI(Tl) crystal is optically coupled to two photomultiplier tubes (PMTs) with the signals recorded as 8-μs waveforms. A trigger is generated when a signal corresponding to one or more photoelectrons occurs in each PMT within a 200-ns time window ().

The analysis presented here uses 1.7 years of data, previously used for the first annual modulation search (), and background modeling with 1-keVee energy threshold (). The data were acquired between 21 October 2016 and 18 July 2018. Three of the eight crystals were observed to have high noise rates in the region of interest (ROI) and so were excluded from the analysis, resulting in an effective data exposure of 97.7 kg × year (, ).

In the ROI, it was found that PMT noise causes most of the triggers. A multivariable boosted decision tree (BDT) () was used to characterize the pulse shapes to discriminate these PMT-induced noise events from radiation-induced scintillation events (, ). To improve the discrimination power, a likelihood score was introduced as an input training variable to the BDT, which rates how well the waveform matches either scintillation events or PMT-induced noise events. The likelihood score particularly enhances the removal of noise pulses and allows us to achieve a 1-keVee threshold (). The BDT is trained with samples of scintillation-rich 60Co calibration data and PMT noise–dominant single-hit physics data. The multiple-hit events consist of hits in multiple crystals or liquid scintillator that cannot be caused by dark matter interactions. The event selection efficiencies for scintillation events are evaluated with the 60Co calibration dataset and cross-checked with the physics data, as well as nuclear recoil events. The efficiencies from the 60Co calibration data were found to be consistent with previously measured efficiencies for nuclear recoil events obtained using a monoenergetic 2.42-MeV neutron beam () as shown in Fig. 1. The efficiency differences and their uncertainties are included as a systematic uncertainty.

Fig. 1.

Efficiencies for β/γ and nuclear-recoil events.

Blue dots show the efficiencies for β/γ events for a COSINE-100 crystal. Black and red dots are efficiencies of β/γ and nuclear recoil events, respectively, for a small-size test crystal. This test crystal was cut from the same ingot of the COSINE-100 crystal and used for the neutron beam measurement. All measurements are consistent within the systematic uncertainty of the efficiency shown in the gray band.

Background modeling

Events in the remaining dark matter search dataset predominantly originate from environmental γ and β radiations. Sources include radioactive contaminants internal to the crystals or on their surfaces, external detector components, and cosmogenic activation (). To understand this, the background spectrum for each crystal is modeled using computer simulations based on the Geant4 toolkit ().

Events are classified according to their energy: 1 to 70 keVee are low energy and 70 to 3000 keVee are high energy. The single-hit and multiple-hit data are separated in the background modeling of the NaI(Tl) crystals. To understand background spectra, Geant4-based simulation events are generated and recorded in a format that matches that of the COSINE-100 data acquisition system. Energy resolutions and selection efficiencies for each crystal are applied. The fraction of each background component is determined from a simultaneous fit to the four measured distributions. For the single-hit events, only 6 to 3000 keVee events are used to avoid a bias of the WIMP signal in the ROI. Details of the background modeling for the dataset are described elsewhere ().

The background components are divided into four categories: internal contamination, surface contamination, external sources, and cosmogenic activation. The 238U, 232Th, 40K, and 210Pb contaminations in the crystal constitute the internal background. The 210Pb contaminations on the crystal surface and adjacent materials are the surface component. Backgrounds from 238U, 232Th, and 40K in the PMTs, liquid scintillator, and the shield materials constitute external sources. To estimate contributions from cosmogenic activation, we use a time-dependent analysis that takes into account the cosmic-ray exposure time on the ground and the cooling time in the underground laboratory of each crystal ().

The most dominant background components in the ROI are generated by internal radionuclide contamination and by cosmogenic activation. This includes 210Pb and 40K internal contaminants and 210Pb surface contamination. The contribution to the ROI from cosmogenic activation is mostly due to 3H with some additional contributions from 113Sn and 109Cd. Background modeling was performed independently for each crystal, and Fig. 2 shows the accumulated result of the model fit to data and the systematic uncertainties.

Fig. 2.

Energy spectra.

Single-hit (A) and multiple-hit (B) data presented here are summed energy spectra for the five crystals (black dots) and their background modeling (red solid line) with the 68 and 95% confidence intervals. The expected contributions to the background from internal radionuclide contaminations, the surface of the crystals and nearby materials, cosmogenic activation, and external backgrounds are indicated. The 1- to 6-keV region of the single-hit spectrum is masked because these events are not used for the background modeling.

Several sources of systematic uncertainties are identified and included in this analysis. The largest systematic uncertainties are those associated with the efficiencies, which include statistical errors in the efficiency determinations with the 60Co calibration and systematic errors derived from the independent cross-checks of the physics data and the nuclear recoil events. Uncertainties in the energy resolution and nonlinear responses of the NaI(Tl) crystals () affect the shapes of the background and signal spectra. The depth profiles of 210Pb on the surface of the NaI(Tl) crystals, studied with a 222Rn-contaminated crystal, are varied within their uncertainty (). Variations in the levels and the positions of external uranium and thorium decay-chain contaminants are also taken into account. Effects of event rate variations and possible distortions in the shapes of spectra are considered in systematic uncertainties.

Dark matter interpretations

We consider various WIMP models to determine the possible contribution from WIMP interactions to the measured energy spectra using the simulated data. The DAMA/LIBRA-phase2 data () were found not to be compatible with the canonical model (, ), which is an isospin-conserving spin-independent interaction between WIMP and nucleus in the specific context of the standard WIMP galactic halo model and is the most commonly used interpretation of the direct detection of the WIMP dark matter (). However, an isospin-violating interaction in which the WIMP-proton coupling is different from the WIMP-neutron coupling provides a good fit to the observed annual modulation signals from the DAMA/LIBRA-phase2 data (, ). To interpret the DAMA/LIBRA data and compare with the COSINE-100 data, we use the best-fit values of the effective coupling of WIMPs to neutrons and to protons (fn/fp) obtained for the simultaneous fit of DAMA/LIBRA-phase1 and DAMA/LIBRA-phase2 data described elsewhere (). We also interpret the results of the COSINE-100 data in the canonical model for the comparison with the DAMA/LIBRA-phase1–only data.

We use the nuclear recoil quenching factor (QF) from recent measurements with monoenergetic neutron beams () (QF is the ratio of the scintillation light yield from sodium or iodine recoil relative to that for electron recoil for the same energy). In those measurements, neutron-tagging detectors at a fixed angle relative to the incoming neutron beam direction provide unambiguous knowledge of deposited energy. We obtained a strong energy dependence of the nuclear recoil QFs. Modelings of the QF measurements described in () are appropriated for this analysis (subsequently referred to as new QF). However, most studies interpreting the DAMA/LIBRA’s results have used substantial larger QF values that were reported by the DAMA group in 1996 () (subsequently referred to as DAMA QF), which were obtained by measuring the response of NaI(Tl) crystals to nuclear recoils induced by neutrons from a 252Cf source. The best description of the measured nuclear recoil spectra from the 252Cf source obtained 30 ± 1% and 9 ± 1% for the sodium and iodine QF values, respectively, assuming no energy dependence of the QF values. Those values obtained from the new measurements are approximately 13 and 5% for the sodium and iodine, respectively, at 20 keVnr, where keVnr is nuclear recoil energy (). Efficient noise rejection and correct evaluation of trigger and selection efficiencies are essential for proper estimation of the QFs (, , ). Although the measurements of the DAMA QF values were required to check the efficiency evaluations and no energy-dependent QF assumption (), the hypothesis of different QFs () in the NaI(Tl) crystals used by DAMA/LIBRA and COSINE-100 needs to be checked. Note that results from analysis of the previous 59.5 days of COSINE-100 data with a 2-keVee threshold were not sufficient to exclude all the DAMA/LIBRA 3σ regions when different QFs are used ().

To search for evidence of a WIMP signal in the data, we use a Bayesian approach with a likelihood function based on Poisson probability. The likelihood fit is performed to the measured single-hit energy spectra between 1 and 15 keVee for each WIMP model for several masses. Each crystal is fitted with a crystal-specific background model and a crystal-correlated WIMP signal for the combined fit by multiplying the five crystals’ likelihoods. Means and uncertainties for background components, which are determined from the modeling (), are used to set Gaussian priors for the background. The systematic uncertainties are included in the fit as nuisance parameters with Gaussian prior (see Materials and Methods).

A good fit to the DAMA/LIBRA-phase2 data was obtained with the isospin-violating interaction (, ). We simultaneously use the DAMA/LIBRA-phase1 and phase2 data to fit three parameters: the WIMP mass, the WIMP-proton cross section, and fn/fp. The best fits were obtained for two different values of fn/fp, favoring WIMP-sodium and WIMP-iodine interactions as fn/fp = −0.76 and −0.71, respectively. For the best-fit values of fn/fp, the 3σ allowed regions in the WIMP-mass and the WIMP-proton cross-sectional parameter spaces are obtained ().

The COSINE-100 data are fitted to each of the different WIMP masses for each fn/fp value using only the new QF values. An example of a maximum likelihood fit for a 11.5 GeV/c2 WIMP and fn/fp = −0.76 WIMP signal is presented in Fig. 3. The summed event spectrum for the five crystals is shown together with the best-fit result. For comparison, the expected signal for a 11.5 GeV/c2 WIMP with a spin-independent WIMP-proton cross section of 2.5 × 10−2 picobarn (pb), the central value of the DAMA/LIBRA best fit using the DAMA QF values for the WIMP-sodium interaction, is shown by the red solid line. No excess of events that could be attributed to WIMP interactions is found for the considered WIMP signals. The posterior probabilities of signals are consistent with zero in all cases, and 90% confidence level limits are determined (see fig. S4). Figure 4 shows the 3σ contours of the DAMA/LIBRA data in the best-fit values of fn/fp using the new QF values and the DAMA QF values together with the 90% confidence level upper limits from the COSINE-100 data using the same fn/fp and the new QF values. The 90% confidence level limits from the 1.7-years COSINE-100 data show approximately one-order of magnitude better limits than those of our previous results using 59.5-days data and exclude the DAMA/LIBRA allowed 3σ regions for both of the two different QF values.

Fig. 3.

Example fit results for a 11.5 GeV/c2 WIMP mass in the case of fn/fp=−0.76.

Presented here is the summed energy spectrum for the five crystals (black filled circles shown with 68% confidence level error bars) and the best fit (blue line) for which no WIMP signals are obtained. Fitted contributions to the background from internal radionuclide contaminations, the surface of the crystals and nearby materials, cosmogenic activation, and external backgrounds are indicated. The green (yellow) bands are the 68% (95%) confidence level intervals of the systematic uncertainty obtained from the likelihood fit. For presentation purposes, we indicate the signal shape (red line) assuming a WIMP-proton cross section of 2.5 × 10−2 pb corresponding to the DAMA best-fit value for the WIMP-sodium interaction using the DAMA QF values.

Fig. 4.

Exclusion limits on the WIMP-proton spin-independent cross section for the isospin-violating interaction.

The 3σ allowed regions of the WIMP mass and the WIMP-proton cross section associated with the DAMA/LIBRA-phase1 + phase2 data (blue solid contours) using the new QF values in their best fit for (A) sodium scattering and (B) iodine scattering hypotheses are compared with the 90% confidence level exclusion limits from the COSINE-100 data (black solid line), together with the 68 and 95% probability bands for the expected 90% confidence level limit assuming the background-only hypothesis. The dashed blue contours show the allowed regions of the DAMA/LIBRA-phase1 + phase2 data using the DAMA QF values. For comparison, limits from the initial 59.5-day COSINE-100 data () are shown by the purple solid line. In each plot, we fix the effective coupling ratios to neutrons and protons fn/fp to the best-fit values of the DAMA data.

Although the DAMA/LIBRA-phase2 data do not fit with the canonical model, their phase1 data have been shown to be well fit with an isospin-conserving spin-independent WIMP-nucleon interaction (, ). The 90% confidence level upper limits from the COSINE-100 data for the canonical model are also obtained. Figure 5 shows the 3σ allowed regions that are associated with the DAMA/LIBRA-phase1 signal using the new QF values and the DAMA QF values, together with the 90% confidence level upper limits from the COSINE-100 data using the new QF values. These limits mostly exclude the DAMA/LIBRA allowed region although different QF values are considered for each experiment.

Fig. 5.

Exclusion limits on the WIMP-nucleon spin-independent cross section of the isospin-conserving interaction.

The observed (filled circles with black solid line) 90% confidence level exclusion limits on the WIMP-nucleon spin-independent cross section from the COSINE-100 are shown together with the 68 and 95% probability bands for the expected 90% confidence level limit, assuming the background-only hypothesis. The limits are compared with a WIMP interpretation of the DAMA/LIBRA-phase1 3σ allowed region using the new QF (blue solid contours) and the DAMA QF (blue dashed contours) ().

In addition, we have checked each operator in nonrelativistic effective field theory models where previous null results from the 59.5-day COSINE-100 data do not fully cover the 3σ regions of the DAMA/LIBRA data for a few operators (). The 1.7-year data are now found to fully cover the 3σ allowed regions assuming same DAMA QF values as one can be seen in Fig. 6.

Fig. 6.

Exclusion limits on the WIMP-proton cross section for the effective field theory operators.

DAMA/LIBRA 3σ allowed regions (blue contours) and COSINE-100 90% confidence level exclusion limits of previous analysis (pink solid lines) and this work (black dots and lines) on the WIMP-proton cross sections for the effective field theory operators using same DAMA QF values are presented. For each operator, fn/fp is fixed to the corresponding best-fit value of the DAMA/LIBRA data.

DISCUSSION

After the release of the initial 59.5-day COSINE-100 data with null observations using the same NaI(Tl) target material, a few possibilities were suggested that preserve the consistency between the DAMA/LIBRA and COSINE-100 results (, , ). The results of this analysis, with 1.7-year-accumulated COSINE-100 data and improved analysis technique with 1-keVee energy threshold, do not favor these suggested possibilities. A model-independent data analysis of the annual modulation with several-year COSINE-100 data is required for an unambiguous conclusion; nevertheless, these results provide strong constraints on the dark matter interpretation of the DAMA/LIBRA annual modulation signals with the same NaI(Tl) target materials.

MATERIALS AND METHODS

Event selection

An event satisfying the trigger condition of coincident photoelectrons in both of the crystal’s readout PMTs within 200 ns is acquired with 500-MHz flash analog-to-digital converters and is recorded as an 8-μs-long waveform starting 2.4 μs before occurrence of the trigger (). In the offline analysis, muon-induced events are rejected when the crystal hit event occurs within 30 ms after a muon candidate event in the muon detector (, ). In addition, we require that leading edges of the trigger pulses start later than 2.0 μs after the start of the recording, waveforms from the crystal contain more than two single photoelectrons, and the integral waveform area below the baseline does not exceed a limit. These criteria reject muon-induced phosphor events and electronic interference events. A multiple-hit event has accompanying crystal signals with more than four photoelectrons in an 8-μs time window or has a liquid scintillator signal above an 80-keVee threshold within 4 μs (). A single-hit event is classified as one where the other detectors do not meet these criteria.

During the 1.7-year data-taking period, no significant environmental abnormality or unstable detector performance was observed. The high light yield of the six crystals, approximately 15 photoelectrons/keVee, allowed an analysis threshold of 2 keVee in the previous analysis. However, the other two crystals had lower light yields and required higher analysis thresholds of 4 and 8 keVee, respectively (, ). Because their direct impact on the WIMP search is not substantial, we do not include single-hit events from these two crystals in the WIMP search analysis, although they were used in the identification of multiple hits.

In the low-energy signal region below 10 keVee, PMT-induced noise events predominantly contribute to the single-hit physics data in two different ways. The first class has a fast decay time of less than 50 ns compared with a typical NaI(Tl) scintillation of about 250 ns. The second class, which occurs less often than the first, has different characteristics of slow rise and decay time, as characterized in (, ). Noise events of the second class are intermittently produced by certain PMTs. We have developed monitoring tools for data quality verification, including monitoring event rates of the second class of noise. If a crystal had an increased rate due to the second class of noise, then the relevant period of data was removed. One crystal detector had this class of noise for the whole period, but for the other five detectors, more than 95% of the recorded data could be used without the second-class noise-induced events.

A BDT was developed to separate scintillation signals from the first class of noise. The fast PMT-induced events with energies greater than 2 keVee were efficiently removed by the BDT, which is based on multiple parameters (, ): the balance of the deposited charge from two PMTs, the ratio of the leading-edge (0 to 50 ns) to trailing-edge (100 to 600 ns) charge, and the amplitude-weighted average time of the signal. However, using this BDT, the scintillation events with energies below 2 keVee were contaminated by an exponential increase in noise events.

We developed new parameters to characterize the PMT-induced noise events for efficient selection of the scintillation events below 2 keVee. Two likelihood parameters for an event are constructed by templates derived from data samples enriched alternatively in scintillation-signal events and noise-signal events. A 60Co source that produces low-energy electron signals through Compton scattering is used to generate the signal-enriched sample. Fast-decay events in the single-hit data sample are used as the noise-enriched sample.

The likelihood parameter of the event is calculated aswhere T and W are the heights of the ith time bin in the reference template and event, respectively. Likelihood parameters for the scintillation signal events (lnℒs) and the PMT-induced noise events (lnℒn) are evaluated for each event. With these, we define a likelihood score aswhere large p for an event implies that the event is more closely matched to the scintillation signal than the PMT-induced noise events. The updated BDT is trained with the likelihood score. This provided good discrimination against PMT-induced noise events and enabled us to lower the threshold to 1 keVee (). A BDT score, a single discriminating parameter created by combining the various selections for input parameters according to their corresponding importance, for the single-hit physics data near the energy threshold (1 to 1.25 keVee) is presented in fig. S1. With the established selection criteria, we reduced the PMT-induced noise contamination level to less than 0.5%.

Event selection efficiencies for the electron recoil events were evaluated with the 60Co calibration data. The efficiency is calculated by integrating the model distribution for the scintillation signals and the PMT-induced noise events shown in fig. S1. A specialized apparatus with a monoenergetic 2.42-MeV neutron beam was used to measure the selection efficiencies of the nuclear recoil events (). This measurement was performed with a small-size test crystal that was cut from the same ingot as a crystal used for the COSINE-100 experiment. The efficiencies determined for the electron recoil events and the nuclear recoil events are consistent within the 5% level as shown in Fig. 1. Systematic uncertainties in the efficiency measurements account for deviations from different measurements and from the 60Co calibration data and the single-hit physics data.

Systematic uncertainties

In addition to the statistical uncertainties in the background and signal models, various sources of systematic uncertainties are taken into account. Errors in the selection efficiency, the energy resolution, the energy scale, and background modeling technique translate into uncertainties in the shapes of the signal and background probability density functions, as well as to rate changes. These quantities are allowed to vary within their uncertainties as nuisance parameters in the likelihood fit.

The most influential systematic uncertainty is the error associated with the efficiencies shown as the shaded region in Fig. 1. This is because the efficiency systematic uncertainty maximally covers the statistical uncertainties in the ROI. We include two types of systematic variations in the efficiency error. At first, we account maximum variations of the uncertainty across the width of the energy bin to provide ±1σ ranges as shown in fig. S2A. Relative uncertainties shown in fig. S2B are included as nuisance parameters of the background and signal in the likelihood function (see Eqs. 8 and 9). However, because of the dominant errors in the ROI, we consider maximum shape distortions that can mimic the WIMP signal as shown in fig. S2C. Its relative uncertainty is shown in fig. S2D and included as a nuisance parameter.

In the background model fit, background activities are constrained by Gaussian constraint terms added to the likelihood function as determined by measured activities and their uncertainties. The systematic uncertainties associated with the background modeling include the uncertainties of the activities estimated by the background model fit. In addition, different locations of external radioactive contaminations are taken into account by generating external contributions at different positions. Background contributions from 210Pb contamination on the surface of the NaI(Tl) crystals were studied with a small NaI(Tl) crystal exposed to 222Rn from a 226Ra source (). Depth profiles from two exponential components were modeled to fit the 222Rn-contaminated crystal and matched to the test-setup data (). Uncertainties in the measured depth profiles are propagated into systematic uncertainties. Here, we generate the background model associated with each systematic variation and account relative uncertainty to be added as the nuisance parameter.

The energy calibration is performed by tracking the positions of internal β and γ peaks from radioactive contaminations in the crystals, as well as with external sources (). The nonlinear detector response of the NaI(Tl) crystals () in the low-energy region is modeled with an empirical function across all crystals (). Subtle differences for each crystal from the general nonlinearity model of the NaI(Tl) crystals are evaluated to consider the systematic uncertainty on the energy scale. The energy resolution for each crystal is evaluated during the calibration process. In particular, tagged 0.8- (22Na) and 3.2-keVee (40K) x-ray lines in the multiple-hit data are used to determine the energy resolution for low-energy events. Statistical uncertainties associated with the number of tagged events are regarded as the resolution systematic. Figure 2 shows a comparison of the background model to the data together with ±1σ and ±2σ bands of the systematic uncertainties that are evaluated from the quadrature sum of each systematic component. The relative errors are used for the nuisance parameter of the background.

Expected WIMP signal

The differential nuclear recoil rate per unit target mass for elastic scattering between WIMPs of mass mχ and target nuclei of mass Μ is ()where ρχ is the local dark matter density, Enr is the nuclear recoil energy, σ(Μ, Enr) is the WIMP-nucleus cross section and f(v, t) is the time-dependent WIMP velocity distribution. The reduced mass μ is defined as mχΜ/(mχ + Μ) and the minimum WIMP velocity vmin is .

For the WIMP velocity distribution, we assume the SHM ()where Nesc is a normalization constant, vE is the Earth velocity relative to the WIMP dark matter, and σv is the velocity dispersion. The SHM parameterization is used with the local dark matter density ρχ = 0.3 GeV/cm3, vE = 232 km/s, , and the galactic escape velocity vesc = 544 km/s ().

The effective field theory operators and nuclear form factors described in (–) are used to model the nuclear responses in the differential cross section. The generalized spin-independent response () is used for both the isospin-conserving and the isospin-violating spin-independent (SI) interactions. For isospin-violating SI interactions, the WIMP-nucleon coupling coefficient ratio, fn/fp, is fixed to the best-fit values for the DAMA/LIBRA data (). These nuclear responses, including form factors, are implemented using the publicly available DMDD package () to evaluate the WIMP-nucleus cross section, σ(Μ, Enr) in Eq. 3 (raw signal spectra). We subsequently apply the QFs, energy resolution, and selection efficiency to obtain the expected nuclear recoil rate in electron-equivalent energy for the detector, dR/dEee. Figure S3 shows dR/dEnr (A) and dR/dEee (B) spectra for three different WIMP models assuming a WIMP-proton cross section equal to 1 pb.

In addition, the WIMP signals are generated in the context of the nonrelativistic effective theory of WIMP-nucleus scattering that has been tested using previous data without full coverage of the DAMA/LIBRA 3σ allowed regions (). For simplicity, we assume that one of the effective operators allowed by Galilean invariance dominates in the effective Hamiltonian of a spin-half dark matter particle at a time. We use the best-fit neutron-over-proton coupling ratio of the DAMA/LIBRA-phase2 data assuming the DAMA QF reported in (), for each operator. The DMDD package () is also used to generate signal spectra by the effective operators. A few operators are not evaluated because the DMDD package does not include form factors for these operators. Here, we assume the same DAMA QFs for the COSINE-100 data for a simple comparison.

Bayesian approach

A Bayesian approach is adopted to extract the WIMP signal from the COSINE-100 data. For each WIMP interaction model, a posterior probability density in terms of the WIMP-proton cross section σWIMP is obtained from the Bayes’ theorem using a marginalization of the likelihood function that includes the prior ()where P(σWIMP∣Μ) is a probability density function (PDF) and ℒ(Μ∣σWIMP, α, β) is the likelihood function. The prior π(σWIMP) has a constant and zero values for positive and negative σWIMP, respectively. The normalization constant N makes the integration of the posterior PDF to be unity and the Μ represents the measured data. The α and the β denote the nuisance parameters to control the effect by systematic uncertainties. Because the measurements are independent and follow Poisson probabilities, the likelihood function is built by a product of Poisson probabilitieswhere i and j denote the crystal number and the energy bin, respectively, and M is the number of observed events for crystal i in jth energy bin. The π(α, β) is a term for constraining nuisance parameters by systematic uncertainties.

The expected number of events, denoted as E(σ, α, β), is obtained as the sum of the number of events in the signal and backgroundwhere the number of background events B(α, β) and signal events S(σ, α) are generated from the simulated experiments through the background modeling and the WIMP signal discussed above, with effects by systematic uncertainties. The systematic uncertainty affecting the background model is included as a function of the nuisance parameter α and β, aswhere is the number of background events obtained by modeling. The nuisance parameter α controls the effect of the energy-dependent uncertainty, ε, which is 1σ relative error for kth systematic uncertainty. Meanwhile, another nuisance parameter, β, adjusts the activity for lth background component. Similar impact on the WIMP signal is considered by means of the expressionwhere M and T denote the mass and data exposure for crystal i and R is the expected rate of WIMP-proton interaction through an integration of dR/dEee in the jth energy bin. Each nuisance parameter is constrained with evaluated uncertainty assuming a Gaussian distributionwhere δ is the uncertainty of activity for lth background component. The Markov Chain Monte Carlo (, ) via Metropolis-Hastings algorithm (, ) is used for the multivariable integration in posterior PDF. We developed our own Bayesian tool for this process. A comparison with a publicly available Bayesian analysis toolkit () was done for the initial 59.5-day COSINE-100 data, and both tools showed consistent results.

To avoid biasing the WIMP search, the fitter was tested with simulated event samples. Each piece of experimental data is prepared by Poisson random extraction from the modeled background spectrum (), assuming a background-only hypothesis. Marginalization to obtain the posterior PDF for each simulation sample is performed to set the 90% confidence level exclusion limits as shown in fig. S4. The 1000 simulated experiments result in 68 and 95% bands of the expected limit presented in Figs. 4 and 5. The data fits are done in the same way as the simulated data. Figure S4 shows the posterior PDFs and their cumulative distribution functions (CDFs) of data for two different WIMP models. The CDF provides the 90% confidence level exclusion limit for each fit.