Three-year annual modulation search with COSINE-100

COSINE-100 is a direct detection dark matter experiment that aims to test DAMA/LIBRA's claim of dark matter discovery by searching for a dark matter-induced annual modulation signal with NaI(Tl) detectors. We present new constraints on the annual modulation signal from a dataset with a 2.82 yr livetime utilizing an active mass of 61.3 kg, for a total exposure of 173 kg$\cdot$yr. This new result features an improved event selection that allows for both lowering the energy threshold to 1 keV and a more precise time-dependent background model. In the 1-6 keV and 2-6 keV energy intervals, we observe best-fit values for the modulation amplitude of 0.0067$\pm$0.0042 and 0.0051$\pm$0.0047 counts/(day$\cdot$kg$\cdot$keV), respectively, with a phase fixed at 152.5 days.

COSINE-100 is a direct detection dark matter experiment that aims to test DAMA/LIBRA's claim of dark matter discovery by searching for a dark matter-induced annual modulation signal with NaI(Tl) detectors. We present new constraints on the annual modulation signal from a dataset with a 2.82 yr livetime utilizing an active mass of 61.3 kg for a total exposure of 173 kg·yr. This new result features an improved event selection that allows for both lowering the energy threshold to 1 keV and a more precise time-dependent background model. In the 1-6 keV and 2-6 keV energy intervals, we observe best-fit values for the modulation amplitude of 0.0067±0.0042 and 0.0051±0.0047 counts/(day·kg·keV), respectively, with a phase fixed at 152.5 days.

I. INTRODUCTION
Cosmological observations indicate that more than a quarter of our Universe's mass-energy exists in the form of a massive, non-luminous component known as dark matter [1]. The apparent abundance of dark matter has given rise to several experiments over the past three decades that aim to observe dark matter directly [2][3][4][5][6][7][8][9][10]. Despite this large-scale effort, none of these experiments have seen any signal indicating the existence of dark matter, except for the DAMA experiments, DAMA/NaI [11] and DAMA/LIBRA [12,13].
DAMA's claim of dark matter discovery comes in the form of an annual modulation in the event rate of thallium-doped sodium iodide (NaI(Tl)) detectors, which has persisted for more than two decades. This observed * Corresponding author: william.thompson@yale.edu modulation possesses the characteristics expected of a dark matter-induced annual modulation signal [14,15], including a maximum rate near June 2 nd and a period of one year. Additionally, this modulation signal is observed at an extremely high significance of 12.9σ in the 2-6 keV energy range and 9.5σ in the 1-6 keV energy range compared with the no-modulation hypothesis [13]. Despite this high significance, the DAMA result is in severe tension with the null results obtained by other direct-detection dark matter experiments within most commonly considered models of dark matter [16].
The COSINE-100 collaboration aims to resolve this tension by performing a model-independent test of the DAMA collaboration's claim of dark matter discovery [17]. This model independence is achieved by using the same target material as the DAMA experiments, NaI(Tl). Previously, we have excluded spin-independent WIMP-nucleus interactions as the origin of DAMA's modulation signal within the context of the standard halo model, becoming the first NaI(Tl)-based experiment to test the DAMA result [18,19]. We have also published results on a model-independent search for the DAMAobserved modulation with the initial 1.7 yr of data from COSINE-100 [20]. This first modulation search was statistically limited and thus compatible with both the no annual modulation hypothesis and the DAMA experiments' best-fit modulation amplitude in the 2-6 keV region. In addition, the ANAIS-112 collaboration, which also aims to test DAMA's discovery claim with an array of NaI(Tl) detectors, has released results based on one and a half [21] and three years of detector operation [22], the latter of which is incompatible with the DAMA/LIBRA result at 3.3σ.
In this paper, we present the results of a search for a dark matter-induced annual modulation signal in COSINE-100 from October 21, 2016 to November 21, 2019, with a total livetime of 2.82 yr. In addition to the increased exposure, this analysis features a number of improvements to our previous modulation search, including a decrease in energy threshold to 1 keV [23], a timedependent background model based on our modeling reported in [24], and the implementation of a Bayesian analysis approach.
We first provide, in Section II, an overview of the COSINE-100 detector and its performance over the three-year period of data-taking investigated in this analysis. Section III details the updated event selection used to lower our energy threshold to 1 keV. In Section IV, we describe the background model used in this analysis. In Section V, we present our modulation search procedure and results from this analysis, with concluding remarks given in Section VI.

A. Detector Design
A full description of the COSINE-100 detector has been previously published in Refs. [17,25]. In this section, we provide a brief overview of aspects of the detector configuration and performance that are relevant to the annual modulation search.
The COSINE-100 detector consists of eight lowbackground NaI(Tl) detectors, with a total mass of 106 kg, located at a depth of approximately 1800 meters water equivalent in the Yangyang underground lab in South Korea. These detectors are submerged in a 2200 L liquid scintillator (LS) detector that serves as a background veto system [26]. This veto reduces backgrounds from radiation external to the NaI(Tl) detectors, in addition to enabling the tagging and removal of internal low-energy backgrounds that are accompanied by a highenergy emission, particularly decays from 40 K. The LS is continually flushed with nitrogen gas in order to remove radon from the detector volume and prevent discoloring of the scintillator from contact with oxygen. External backgrounds are further reduced by a shielding structure composed of a layer of 3-cm thick copper surrounded by 20 cm of lead. Thirty-seven 3-cm thick plastic scintillator panels placed around the shielding structure provide 4π cosmic ray muon tagging [27,28].
Each NaI(Tl) crystal detector is optically coupled to two photomultiplier tubes (PMTs) that detect scintillation photons from energy depositions in the detector. These eight detectors are referred to as Crystal 1 (C1) through Crystal 8 (C8). C1, C5, and C8 are excluded from this analysis due to high background levels and low light yields in the case of C5 and C8, and a high noise rate in the case of C1. This results in a total effective mass of 61.3 kg.

B. Detector Stability
The COSINE-100 detector room is instrumented with sensors that monitor the environmental conditions of the room and the detector. This includes sensors that monitor the room temperature, humidity, and radon level, and the temperature of the LS volume, which serves as a proxy for the NaI(Tl) detectors' temperatures. Continuous monitoring of these environmental conditions allows us to assess the stability of the environment and detector in real-time. The monitoring system also allows us to perform offline investigations into possible correlations between environmental conditions and the detector event rate. A detailed description of the monitoring system can be found in Ref. [29].
We calibrate the NaI(Tl) detectors using both external γ-ray sources and β-and γ-ray emissions from radioactive contaminants internal to the NaI(Tl) detectors. At energies below 70 keV, the nature of the non-linear light response of NaI(Tl) is studied using low-energy internal decays. We account for this non-linearity in our energy calibration via an empirical model based on these low-energy decays. Full details of the energy calibration method can be found in Ref. [24].
In addition to accounting for non-linearities in the light yield at low energies, it is necessary to monitor and correct for gain shifts in the NaI(Tl)-coupled PMTs to ensure a stable detector energy scale. Stability of the energy scale is of particular importance for annual modulation searches, as instabilities can manifest as changes in the event rate, possibly inducing an artificial annual modulation signal. To guard against this effect, we monitor and correct for gain shifts in each of the NaI(Tl)-coupled PMTs by tracking changes in the reconstructed energy of the ∼50 keV decay from 210 Pb. This decay line originates from the combination of an electron emitted in the β-decay of 210 Pb to 210 Bi (Q = 17 keV) and a prompt 46.5 keV γ-ray from the corresponding relaxation of the excited 210 Bi nucleus. Since this decay originates from radioisotopes internal to the NaI(Tl) detector, it enables continuous monitoring of the PMT gains over the course of the data run. In practice, we measure the reconstructed energy of this decay peak over contiguous 20-   hour periods over the full course of the run. The average size of the gain shift across the NaI(Tl)-coupled PMTs is < 10% over the three-year-long dataset and is less than 15% in all PMTs. To stabilize each NaI(Tl) detector, we apply a time-dependent correction factor to the energy scale of each PMT based on the observed evolution of the reconstructed energy of the 210 Pb peak over time.
After the application of this correction, we assess the stability of each gain-corrected NaI(Tl) detector in the 1-6 keV region of interest (ROI) by tracking the measured energy of the 3.2 keV x-ray from the decay of 40 K within each sodium iodide detector over time. As this 3.2 keV x-ray is emitted in coincidence with a higherenergy 1461 keV γ-ray, we are able to select a high-purity sample of 3.2 keV energy depositions in each NaI(Tl) detector by requiring a coincident energy deposition between 900 and 1550 keV within our LS veto. With the gain-shift correction applied, we find that the energy scale over time is stable within the ROI in all five NaI(Tl) detectors used in this analysis, as illustrated in Fig. 1. Specifically, we find no evidence of a change in any of the detectors' energy scales over time, with a maximum reduced chi-square statistic of χ 2 /DOF= 20.2/19 across all five detectors.

III. EVENT SELECTION
In this study, we employ a newly developed event selection procedure that enables lowering the analysis energy threshold to 1 keV, an improvement compared with the 2 keV threshold in our previous modulation search [20]. A full description of the updated event selection procedure can be found in Ref. [23]. Here, we provide a summary of this procedure.
The first step in the event selection procedure is to identify and remove muon-induced events, defined as any event that occurs within 30 ms of an energy deposition in the muon veto system. We also remove events with characteristics that are consistent with noise signals that originate from electronic pickup. Next, we classify the remaining events as "single-hit" or "multi-hit." Specifically, a single-hit event is any event in which there is an energy deposition in only a single NaI(Tl) detector along with no measurable energy deposition in the LS veto; a multi-hit event is any event with an energy deposition in multiple NaI(Tl) detectors, or an event with an energy deposition in at least one NaI(Tl) detector and the LS veto. Dark matter, due to its extremely small probability of interacting with ordinary matter, is expected to exclusively induce single-hit events. While multi-hit events are not directly used in the dark matter search, they provide a convenient sideband sample for the evaluation of cut efficiencies and the assessment of other properties of our modulation search procedure.
The primary focus of our event-selection procedure is to remove PMT-induced noise events that mimic the shape of NaI(Tl) scintillation signals. At lower energies, it becomes increasingly difficult to distinguish between signal and PMT-induced noise events due to the fewer number of photons emitted by a NaI(Tl) detector. Thus, decreasing our energy threshold required the development of a more powerful noise discrimination algorithm than that of our previous approach [18]. This updated event selection procedure utilizes a new metric, named the "likelihood score," that compares individual PMT waveforms to template signal and noise waveforms. The signal template is obtained from a high-purity calibration dataset of scintillation events. The events in this highpurity dataset originate from Compton-scattered γ-rays from a 60 Co calibration source and were isolated from non-scintillation noise events by selecting only multi-hit events. The noise template is created using events identified as noise via an independent pulse-shape identification procedure. The difference between the signal and noise likelihood values is defined as the likelihood score for each individual event.
We utilize a boosted decision tree (BDT) to combine the likelihood score with the event energy and other event discrimination features that were used in our previous modulation search [20]. These other features are used to further quantify PMT waveforms shapes, as well as the asymmetry of light collection between the two PMTs coupled to a single NaI(Tl) detector. Though the precise weighting of each feature by the BDT varies between the five detectors, their general order of importance is the same. The likelihood parameter is assigned the largest weight by the BDT, followed by the event energy. The asymmetry feature is the third highest weighted, with high asymmetry values being characteristic of noise events. The asymmetry feature is followed in importance by parameters characterizing the prompt and delayed charge measured by each of the two PMTs attached to a given NaI(Tl) detector. The other waveform shape parameters used in previous analyses by COSINE-100 make up the remaining, lower-ranked parameters [18,20]. The high ranking of the likelihood parameter illustrates its importance in the updated event selection procedure used in this analysis.
To compute the efficiency of the event selection as a function of energy, we again utilize a high-purity dataset obtained from a 60 Co calibration run. The selection efficiency is defined as the ratio between the number of Compton-scatter events after to the number of Comptonscatter events before the selection and is computed independently for each of the five detectors used in this analysis. We find that this updated event selection achieves a signal selection efficiency above 80% in the 1-1.5 keV energy bin and that approaches unity above 2 keV [19]. For this analysis, the selection efficiency uncertainty of each NaI(Tl) detector was obtained from the statistical uncertainties of the 60 Co dataset. Because of the known difference in waveform shape between electronic and nuclear recoil events in NaI(Tl) detectors [30], we also measure the efficiency of the signal selection on lowenergy nuclear recoil events obtained by irradiating a small NaI(Tl) detector with 2.42 MeV neutrons from a D-D neutron generator [31]. This cross-check is of particular importance given that dark matter is expected to interact with NaI(Tl) detectors by inducing nuclear recoils. The selection efficiencies measured with the nuclear recoil dataset are consistent with those measured using the 60 Co dataset [19], verifying that the event selection does not preferentially remove events generated by dark matter interactions.
Lastly, two-hour periods of data are occasionally removed from consideration in our analysis due to environmental or detector instability. Any large variation in the environmental or detector parameters recorded results in the removal of the corresponding two-hour-long sub-run. One example of such an instability is a discrete spike in the detector trigger rate due to the passage of a high-energy cosmic ray muon through the detector; in this case the corresponding sub-run would be removed to mitigate any possible impact of long-lived muon-induced phosphorescence. In addition, there were a few occasions on which a NaI(Tl) detector experienced a large increase in the number of PMT-induced noise events, resulting in significantly elevated event rates over a sustained period of time [19,20]. Because these periods of high event rates are not correlated between NaI(Tl) detectors, exclusion of a particular sub-run from the analysis occurs on a detector-by-detector basis. While we expect the event selection procedure to remove the large majority of these events, we nonetheless have fully removed these periods of instability from consideration in this analysis to preclude the possibility of observing a noise-induced modulation signal. The main periods affected by these elevated noise rates are from October 21 to December 19, 2016 and July 27 to August 8, 2019 for C2, and January 12 to March 31, 2017 for C7. After accounting for the removal of these sub-runs and detector calibration periods, we achieve a 91.5% livetime with the COSINE-100 detector from the operating period spanning from October 21, 2016 to November 21, 2019, resulting in a total data exposure of 173 kg·yr. Additionally, 5% of the time over this period is composed of deadtime, which is defined as both periods during which the COSINE-100 detector was not recording data and periods excluded due to detector instability. Lastly, detector calibrations comprise 3.5% of the detector operating period considered in this analysis.

IV. BACKGROUND MODEL
As the search for an annual modulation signal is primarily focused on the COSINE-100 detector's timedependent event rate, it is vital to understand any timedependent background components that can affect the signal rate in our 1-6 keV ROI. Mismodeling these timedependent background components can bias estimates of the amplitude of the sought annual modulation signal, as has been discussed in Buttazzo et al. [32] and Messina et al [33].
A primary focus of the analysis reported here is the development of a background model that fully accounts for each of these time-dependent components, improving upon the background model used in our previous modulation search [20]. Primarily, these short-lived components are cosmogenic radioisotopes created by the interaction of sodium and iodine nuclei with high-energy cosmic rays before the NaI(Tl) detectors were brought underground. We identified particular cosmogenic isotopes expected to contribute to the time-dependent background of the experiment based on prior studies of cosmogenically activated NaI(Tl) detectors [34][35][36][37]. We have previously reported an independent measurement of these cosmogenic isotopes' activities [38] and incorporated this information into the full, time-integrated background model of COSINE-100, as detailed in Ref. [24]. Additionally, we have improved the COSINE-100 background model by refining the spatial distribution of 210 Pb contamination in our NaI(Tl) detectors, also detailed in Ref. [24]. As 210 Pb is the primary background component in the majority of the NaI(Tl) detectors, it is an isotope of particular importance to our annual modulation search.
COSINE-100's time-integrated background model was developed by first simulating the predicted energy spectrum from each radioisotope expected to contaminate the detector using Geant4 [39]. The activity of each iso-  2. Best-fit background model compared with observed data (black histogram) averaged across all five NaI(Tl) detectors used in this analysis from Ref. [24]. The red histogram is the total energy spectrum from our background model. Isotopes with activities modeled as constant in time in our analysis are grouped into the "flat" component (blue histogram). The energy spectra of isotopes whose time dependence is considered within our annual modulation search model are shown individually. The gray shaded region denotes the region of interest for the annual modulation search. tope was then determined by jointly fitting the amplitude of the simulated energy spectra to the energy spectrum measured by a given NaI(Tl) detector. These fits were performed separately for each detector, resulting in independent background models for each detector. Full details of the time-integrated background model are described in Ref. [24]. Fig. 2 displays the total time-integrated background model averaged across the five detectors used in this analysis in "daily rate units," or dru, a common shorthand for counts/(day·kg·keV). The energy spectra of specific isotopes of interest, defined below, are also shown. Worthy of particular note is 3 H, which is the most prominent background component in the ROI in two of the NaI(Tl) detectors used in this analysis; 210 Pb is the most prominent background in the other three detectors. Lshell emissions from cosmogenic isotopes also contribute significantly to the detector background as seen in the figure. The remaining isotopes, whose activities do not vary appreciably on a three-year timescale, are grouped together in a constant-in-time, or "flat," component. The flat component consists primarily of decays of 40 K and isotopes belonging to the 238 U and 232 Th decay chains.
The time-dependent background model is derived directly from this time-integrated background model. First, we determine the average event rate of each timedependent background component within a given energy range by integrating its simulated spectrum over said energy range. Next, we compute the event rate at the beginning of the dataset considered in this analysis induced by each time-dependent component using its average event rate, half-life, and livetime history. From this calculation, we identify eight radioisotopes that are expected to have a measurable impact on the timedependent event rate of COSINE-100. These are defined as isotopes whose decays induced >0.01% of the total event rate in any of the five NaI(Tl) detectors on October 21, 2016, which corresponds to an event rate of 3 × 10 −4 dru. We also required that these isotopes have short half-lives, conservatively defined as less than one million years. The specific isotopes identified and their average activities over the course of the dataset used in the three-year modulation search are reported in Table I. Seven of these isotopes are of cosmogenic origin, highlighting the importance of their inclusion in our background model. We define the total time-dependent background model of each NaI(Tl) detector as the sum of the constant background component and the eight exponentially decaying background components, with the activity of each time-dependent component on October 21, 2016 treated as a nuisance parameter in the annual modulation search.

A. Modulation Search Procedure
We define the total time-dependent event rate of the i th NaI(Tl) detector as the sum of the event rate predicted by the time-dependent background model added to a dark matter-induced, annually modulating signal component (1) In Eq. 1, the sum over k represents the sum over the N bkgd = 8 background components with decay constants λ k . α i is the rate from the constant background component and β i k is the initial rate from the k th background component, both in the i th detector. The modulation is described by its amplitude S m , phase t 0 , and period T = 365.25 days, where ω = 2π/T . Our model requires the same modulation signal amongst all detectors, but allows for different background activity levels to account for the different contamination levels across different detectors.
In order to efficiently treat the large number of nuisance parameters in our model from these background components we utilize a Bayesian analysis framework, making use of Markov chain Monte Carlo techniques to calculate the marginalized posterior distributions of the parameters of interest. The posterior distribution is given by where N is a normalization constant, π(S m , α, β) represents the prior distributions, and α and β are vectors whose components represent the activities of the constant and time-dependent backgrounds in each crystal, respectively. The observed data is denoted by x and comprises the efficiency-and livetime-corrected event rate in each time bin of each NaI(Tl) detector.
To obtain the efficiency-and livetime-corrected event rate over time of a given detector, we first compute the number of events in each 15-day time bin after application of the event selection procedure described in Sec. III. We then evaluate the detector livetime in each time bin and normalize each bin based on its relative exposure. This process accounts for variations in exposure induced by both detector-off periods and data periods removed due to detector instability. We next scale the event rate in each time bin by the reciprocal of the measured selection efficiency to obtain the sought corrected event rate over time. The uncertainty of each time bin is defined as the square root of the number of events in each time bin before correction, scaled by the efficiency and livetime corrections. Lastly, we found that the uncertainties of the measured selection efficiencies had a negligible impact on the measured modulation amplitudes and phases; they are therefore excluded from the remainder of the analysis.
We utilize a binned likelihood built as the product of Gaussian probabilities where µ ij is the expected number of counts in the j th time bin of the i th detector, obtained by integrating R i (t|S m , α i , β i ) over the duration of the j th time bin, and σ ij is its associated uncertainty. N det = 5 is the number of detectors used in this analysis, and N i bin is the number of time bins in the i th detector.
We assume Gaussian priors for the initial activities of all background components, with the mean and standard deviation of each component's priors set to its computed activity and its associated uncertainty on October 21, 2016, as described in Sec. IV. In this paper, we summarize the marginalized posterior distributions of our parameters of interest, the modulation amplitude and phase, by reporting their highest-density credible intervals (HDIs) for one-dimensional posteriors or highest-density credible regions (HDRs) for two-dimensional posteriors.
In our primary analysis, we fix the phase of the modulation such that the maximum occurs on June 2 nd , 152.5 days from the start of the calendar year, as predicted in the standard halo model [14], and search for an annual modulation signal in the 1-6 keV and 2-6 keV ROIs with a flat prior for the modulation amplitude. We also search for an annual modulation signal by allowing both the amplitude and phase of the signal to vary in the fit, with both assigned flat priors. Lastly, we also perform an analysis of the modulation amplitude as a function of energy by dividing the data into 1-keV wide bins from 1-20 keV and fitting for the modulation amplitude in each bin with the phase of the modulation signal again fixed to 152.5 days; this analysis is also performed on the multi-hit data set. As we expect to observe no significant modulation in either the 6-20 keV single-hit or 1-6 keV multi-hit datasets, these datasets serve as sideband regions to validate our fitting procedure and background model.

B. Pseudo-Experiment Validation
To examine the validity of our fitting procedure, we generate five pseudo-experiment ensembles that are analyzed using the same fitting procedure as is applied to the measured data. For a single pseudo-experiment, the event rate over time in each NaI(Tl) detector is generated from our signal plus background model described by Eq. 1 with Poissonian fluctuations introduced in each 15-day time bin. Each ensemble consists of 1000 pseudoexperiments. The initial activities of each background component are determined from our background model. In order to consider the possible parameter space of initial activities of the background components in our model, we randomly select the initial activities of each component from their respective prior distributions. The amplitude of the injected modulation signal is varied between the five different ensembles to elucidate characteristics of our modulation search procedure across an amplitude range that extends from 0.0000 to 0.0105 dru.
We quantify the validity of our search procedure by investigating the bias distribution of the fitted modulation amplitudes of each pseudo-experiment ensemble. Each bias distribution is fitted with a Gaussian function, from which the mean bias of the ensemble is determined. At the modulation amplitudes investigated, we find a maximum bias in our search procedure of -0.0003 dru, which occurs for a 0.0050 dru input modulation amplitude. The results of this investigation are summarized in Fig. 3. Additionally, we utilize these pseudo-experiment ensembles to evaluate the projected uncertainty of the measured amplitude for modulations of various sizes. This is achieved by investigating the distribution of the uncertainties of the fitted modulation amplitudes in each pseudo-experiment ensemble. Based on this analysis, we expect our uncertainty to fall within the interval (0.0037, 0.0041) dru at a 1σ confidence level; we do not observe any dependence of our projected uncertainty on the amplitude of the modulation. Given that the maximum bias observed in this investigation is an order of magnitude smaller than the projected uncertainty, we do not adjust for this bias in our analysis.
In addition to assessing the bias induced by our full model of the event rate over time, we also investigate the bias induced by the simplified background model used in our previous modulation search. This simplified model consisted of a constant and exponentially decaying component in each detector, in addition to the annual modulation signal. The single exponentially decaying component in each detector served to model short-lived cosmogenic isotopes, similar to the exponentially decaying components in the background model described by Eq. 1; however, the specific cosmogenic components present in the detectors were not known at the time of our previous analysis, precluding the use of the more accurate background model presented here. To quantify the bias in the modulation amplitude measured using this simplified background model on the 173 kg·yr dataset studied here, we fit this simplified model to the aforementioned pseudo-experiments. Across all pseudo-experiment ensembles the magnitude of the mean bias is greater than 0.0085 dru, as shown in Fig. 3. This large bias, with a magnitude roughly as large as the DAMA-observed modulation amplitude, illustrates the importance of developing an accurate time-dependent background model in annual modulation searches.

C. Results and Discussion
With the phase of the modulation signal fixed at 152.5 days, we find a best-fit modulation amplitude in the 1-6 keV energy region of 0.0067±0.0042 dru. The observed event rate over time overlaid with the phase-fixed bestfit model for this energy region is shown in Fig. 4, and the marginalized posterior distribution of the modulation amplitude in the 1-6 keV region is shown in Fig. 5. We also perform a fit to our data in the 2-6 keV region and find a best-fit amplitude of 0.0051 ± 0.0047 dru. The measured modulation amplitudes in both energy regions are consistent with both the DAMA-observed modulation and the case of no observed modulation. Table II provides a full summary of the results from this phasefixed modulation search.
In judging compatibility of a best-fit model with the data, we do not assume Wilks' theorem [40] to compute the p-value but instead numerically compute the chi-square distribution from ensembles of 10000 pseudo-  experiments generated from the best-fit model. As an example, the computed chi-square distribution from the 1-6 keV modulation search and the measured chi-square value of the best-fit model are shown in the right panel of Fig. 6. In this case, we find that the best-fit model agrees with the data at a p-value of p = 0.239. In the 2-6 keV energy region, we find a p-value of p = 0.485. The data is in good agreement with our best-fit models, as evidenced by the calculated p-values. In addition, the level of precision obtained in these measurements is in agreement with predictions of our uncertainty from the pseudo-experiment analysis described in Sec. V B, as can be seen in the left panel of Fig. 6.
We also compare the activities of the background components measured in this study with those measured in our search for spin-independent WIMP-nucleus scattering using 1.7 years of data [19]. We find that the measured activities of both the constant-in-time and shorterlived background components are in good agreement between the two analyses; however, the measured values of the 3 H and 210 Pb activities are in moderate tension. An initial inspection of this effect finds that the 210 Pb present on the surfaces of the NaI(Tl) crystals [24] exhibits a decay time that is longer than expected of 210 Pb, while the 210 Pb contained within the bulk of the crystal does not exhibit this effect. An initial assessment of the impact of this effect on the measured modulation amplitude is performed by fitting the measured event rate over time using the model described by Eq.1 with an additional longer-lived 210 Pb component added to each NaI(Tl) detector. This updated model returns identical best-fit values of the modulation amplitude to those not including this longer-lived 210 Pb component, while also improving the agreement between the 3 H and 210 Pb background components measured in Ref. [19] and this analysis. Further investigations are currently underway to confirm the longer-lived 210 Pb component as the cause of the discrepancy between the two analyses, and to confirm the origin of the extended decay time.
The two-dimensional marginalized posterior distributions obtained from the phase-floated modulation search are shown in Fig. 7  ison. It should be noted that the posterior distributions are biased towards positive modulation amplitudes due to the non-linearity of the phase in the fitted model, as described in Ref. [22]. Nonetheless, we are able to evaluate the level of agreement of the best-fit point for each energy range with the no-modulation case. This is done by comparing the best-fit point of the measured posterior distribution with the posterior distributions of the ensemble of pseudo-experiments with no injected modulation signal. Specifically, we compute the median of the distribution of the maximum amplitude within the 68.3% HDR of the posterior of each pseudo-experiment, with the computed values shown in Fig. 7. Having taken the effect of this bias into account, we find the results from the phase-floated modulation search are in agreement with both the DAMA-preferred modulation signal and the case of no annual modulation, as in the fixedphase search.
Lastly, we present the best-fit modulation amplitude as a function of energy for both single-hit and multihit events in Fig. 8. In the two lowest-energy bins, we find a slight preference for the DAMA-observed modulation amplitude, though the no-modulation case also falls within the 95% HDI of the posterior of both energy bins. This result is consistent with our observations in both the phase-fixed and phase-floated modulation searches in the 1-6 and 2-6 keV energy regions. As described in Sec. V A, no dark matter-induced annual modulation signal is expected in the 6-20 keV singlehit and 1-6 keV multi-hit sideband regions, allowing them to serve as crosschecks of our modulation search procedure. We find that the result in the 6-20 keV single-hit sideband is consistent with no modulation at a level of χ 2 /DOF=12.8/14. Additionally, the multihit 1-6 keV result is consistent with no modulation at χ 2 /DOF= 3.53/5.

VI. CONCLUSIONS
In conclusion, we have performed a search for a dark matter-induced annual modulation signal in NaI(Tl) de-tectors with 2.82 yr of data obtained between October 21, 2016 and November 21, 2019 with COSINE-100. We have improved upon our previous modulation search [20] by implementing a more powerful event selection procedure that allowed us to lower our energy threshold to 1 keV and by implementing a fully featured, timedependent background model based on dedicated background studies of short-lived components of cosmogenic origin in COSINE-100. With the phase and period of the modulation signal fixed, we observe a best-fit modulation amplitude of 0.0067 ± 0.0042 (0.0051 ± 0.0047) dru in the 1-6 (2-6) keV signal region, consistent with both the modulation amplitude reported by DAMA and the no-modulation case. In addition, when allowing both the phase and the amplitude of the modulation signal to float as free parameters we measure a best-fit value that is consistent with both the DAMA-preferred value and the case of no modulation. The validity of our modulation search procedure was confirmed via pseudoexperiment studies and analyses of sideband data samples. Although COSINE-100 is unable to distinguish between the DAMA-observed modulation and no modulation signal after three years of operation, we plan to continue operation of the COSINE-100 detector until at least late 2022, when commissioning for the next phase of the experiment, COSINE-200, is scheduled to commence. Thus, the final exposure of COSINE-100 will increase compared with this analysis by more than a factor of two, significantly improving our sensitivity to DAMA's observed modulation signal.