Inelastic dark matter scattering off Thallium cannot save DAMA

We study the compatibility of the observed DAMA modulation signal with inelastic scattering of dark matter (DM) off of the $0.1\%$ Thallium (Tl) dopant in DAMA. In this work we test whether there exist regions of parameter space where the Tl interpretation gives a good fit to the most recent data from DAMA, and whether these regions are compatible with the latest constraints from other direct detection experiments. Previously, Chang et al. in 2010, had proposed the Tl interpretation of the DAMA data, and more recently (in 2019) the DAMA/LIBRA collaboration found regions in parameter space of Tl inelastic scattering that differ by more than $10\sigma$ from a no modulation hypothesis. We have expanded upon their work by testing whether the regions of parameter space where inelastic DM-Tl scattering gives a good fit to the most recent DAMA data survive the constraints placed by the lack of a DM signal in XENON1T and CRESST-II. In addition, we have tested how these regions change with the main sources of uncertainty: the Tl quenching factor, which has never been measured directly, and the astrophysical uncertainties in the DM distribution. We conclude that inelastic DM scattering off Tl cannot explain the DAMA data in light of null results from other experiments.


Introduction
The nature of dark matter (DM) is one of the longest outstanding problems of modern physics. Today there is evidence of the gravitational effects of DM on length scales ranging from dwarf galaxies to the largest observable structures of the Universe. One of the most promising candidates for DM is weakly interacting massive particles (WIMPs). In addition to feeling the effects of gravity, these particles can interact weakly, which makes it possible to detect their interactions with ordinary matter. Despite the massive experimental efforts that have been made to directly detect signatures of WIMPs recoiling off of nuclei, no conclusive evidence of DM interactions with ordinary matter has been found [3][4][5][6][7][8][9].
However, there is one exception to the null-results of direct detection experiments: For the last 20 years the DAMA collaboration has claimed a positive signal of DM scattering, with an accumulated significance of 12.9σ [10][11][12] after three phases of data-taking. The DAMA/LIBRA (previously, DAMA/NaI) experiment consists of 250 kg of NaI(Tl) scintillators and is located deep underground at the Gran Sasso National Laboratory. Unlike other direct detection searches in which experiments must be virtually background free to detect a few signal events, DAMA takes advantage of the expected time variation of the DM-induced recoil rate in a detector due to the Earth's motion around the Sun [13,14]. The effect of annual modulation in the WIMP recoil rate was first proposed in the 1980's by Drukier, Freese and Spergel [15]. While the likely backgrounds that can mimic the DAMA modulation signal have all been rejected by the collaboration, conventional models for WIMP DM also predict a significant signal from elastic scattering off nuclei in other direct detection experiments, making their null results in strong tension with this interpretation of the DAMA signal, see e.g. Refs. [16][17][18]. DAMA/LIBRA has now published data down to a lower threshold of 1keV [12]. With this new data, Ref. [19] showed that the results of the DAMA/LIBRA experiment are no longer self-consistent under the assumption of vanilla spin independent (isospin conserving) elastic scattering. In addition, recent results from the ANAIS-112 experiment [20], with the same NaI(Tl) target as DAMA, show a signal consistent with no modulation and the modulation amplitude is in ∼ 3σ tension with the signal observed at DAMA.
Thus, it is reasonable to ask if there exist models of DM interactions beyond the standard lore that can explain both the signal in DAMA and the lack of one in all other direct detection experiments. Several alternatives have been proposed, such as DM scattering within a generalized effective field theory framework [21,22], models of mirror DM [23,24] and inelastic DM [25,26]. Inelastic DM (iDM) refers to dark matter that scatters off of target nuclei in direct detection experiments by transitioning to an additional state with a different mass. iDM occurs naturally in several theoretical frameworks, such as certain supersymmetric models [27], magnetic iDM [28], and dark photon mediated DM [29]. A recent study of inelastic DM scattering off Na and I at DAMA explores the sensitivity of future annual modulation experiments to several effective models of DM-nucleon interactions which reduce the tension of the observed DAMA signal with null results from other direct detection experiments [30].
In work from 2010, Chang, Lang, and Weiner [1] proposed that the 0.1% of thallium (Tl) dopant in the DAMA crystals might be responsible for the annual modulation observed in DAMA. The stable Tl isotopes have atomic mass of A = 203 and A = 205, and are thus much heavier than most target nuclei in other direct detection experiments (e.g. for xenon, A = 131). When DM scatters inelastically into a heavier state, there exists parameter space where recoils of a given energy are permitted for the heavier Tl and tungsten (W) targets but kinematically forbidden for the lighter xenon (Xe) target. In particular, the minimum velocity required for DM scattering off of lighter targets (e.g. Xe) can be larger than the maximum DM velocity permitted within the Milky Way halo. Thus, in iDM models where the elastic scattering is suppressed, inelastic scattering off of heavy target nuclei (e.g. Tl) can provide the leading contribution to the DM-nucleus scattering cross-section. In addition to the suppression of scattering off of lighter target nuclei in iDM models, the nuclear form factor for heavy target nuclei such as Tl can provide for the characteristic shape of the nuclear recoil spectrum observed in the DAMA signal. Based on data available at the time, the analysis in Ref. [1] suggested that there exist regions of iDM parameter space which could provide a good fit to the DAMA signal and evade constraints from other direct detection experiments. Also, after the most recent phase of data-taking with higher sensitivity to low-energy events, the DAMA collaboration [2] found regions of parameter space where the expected annual modulation of the signal arising from Tl recoils differs by more than 10σ from the null hypothesis of no modulation. This result is only part of the story, however. A 10σ improvement over the null hypothesis is not necessarily impressive if the model is still not a good fit to the data. It is thus vitally important to ensure that a model comparison to the null hypothesis is also coupled with a goodness-of-fit test.
In this work, we test whether there exist regions of iDM parameter space where the Tl interpretation gives a good fit to the most recent data from DAMA and whether these regions are compatible with constraints from other direct detection experiments, such as XENON1T [5,6] and CRESST-II [8]. We also study how the best-fit signal regions depend on the choice of astrophysical parameters and the Tl quenching factor, both of which are subject to large uncertainties. We test the goodness-of-fit for iDM scattering to the observed modulation signal in two cases: Isospin conserving DM interactions with nucleons and an alternative isospin-violating model where Iodine (I) scattering is maximally suppressed in DAMA. We find that there exist regions in the model parameter space that give a good fit to the DAMA modulation but are ruled out by constraints from XENON1T and/or CRESST-II. Thus, iDM scattering off Tl can no longer reconcile the observed annual modulation in DAMA with the lack of a corresponding signal in other direct detection experiments.
The remaining sections of this paper are organized as follows: In Sec. 2 we give an overview of the calculation of the recoil rate for these dark matter models in the DAMA experiment and discuss the Tl quenching factor. In Sec. 3 we discuss the details of our analysis, particularly the effects of varying the astrophysical parameters on the modulation signal and the constraints from other direct detection experiments. In Sec. 4 we present the best-fit regions of these models to the DAMA signal under a variety of assumptions and in Sec. 5 we discuss our conclusions.

Inelastic Scattering in DAMA
Direct detection experiments search for the signatures of nuclear recoils induced by the interactions of DM with an instrumented volume of target nuclei [31]. When elastic scattering of the target nucleus, X SM , is either forbidden or highly suppressed in models which contain at least two DM particles χ and χ that are nearly degenerate in mass, inelastic scattering processes of the type χ + X SM → χ + X SM can provide the leading contribution to the scattering cross section. Such scenarios feature two DM states, χ and χ , where χ is the incident DM particle and χ is some heavier state. The mass difference between the states is denoted by δ = m χ − m χ .
The minimum speed required for inelastic scattering to yield a nuclear recoil of energy where m X is the mass of the target nucleus and µ is the reduced mass of the DM-nucleus system. In the limit δ → 0, we recover the usual expression for elastic scattering. As we discuss in more detail below, the minimum speed required for DM with mass m χ ∼ 100 GeV to induce nuclear recoils with energies E R δ ∼ O(100) keV can exceed the typical escape velocity for the Milky Way halo, ∼ 550 km/s, for target nuclei significantly lighter than Tl. Thus, the inelastic scattering in direct detection experiments with liquid noble gas targets can be significantly suppressed. Yet, due to the presence of the heavier Tl in the DAMA detector, in combination with effects of the Tl nuclear form factors on the recoil spectrum observed at DAMA, the DAMA signal can be fit well by Tl recoils. The goal is to check whether or not the bounds from detectors made of materials lighter than Tl, albeit weakened by kinematic effects, are still strong enough to rule out the hypothesis of inelastic scattering with Tl as an explanation for the DAMA annual modulation data. The differential recoil rate per unit detector mass for DM-nucleus scattering is where f ( v, t) is the DM velocity distribution function in the detector frame and dσ/dq 2 is the differential scattering cross section with q 2 = 2m X E R as the momentum exchange. We take the local DM density to be ρ χ = 0.3 GeV/cm 3 , although different estimates can vary significantly. 1 For a target consisting of several types of nuclei the total differential recoil rate is given by where X denotes a specific nucleus within the target and ξ X is the associated mass fraction. The differential cross section for the recoil of a given nucleus is given by where σ 0 is the scattering cross section at zero momentum transfer and F is the Helm form factor. The Heaviside step function Θ(q max −q) comes from the maximal momentum transfer allowed for the case of inelastic scattering, given the upper velocity cutoff set by the escape velocity of the DM. In order for a nuclear recoil of a given energy to occur, the maximal DM velocity must exceed the minimum velocity to which the detector is sensitive. Thus the Heaviside function can be replaced by imposing a minimum velocity cutoff, v min , in the integral of Eq. (2.2). For spin-independent scattering, we can parameterize σ 0 as: where Z (A − Z) is the number of protons (neutrons) in the target and f p (f n ) is the coupling of the WIMP to the proton (neutron). In this analysis we consider two cases: Isospin conserving scattering: f n /f p = 1 (2.6a) Isospin violating scattering: f n /f p = −53/(127 − 53) (2.6b) The isospin violating case is tuned so that DM-I scattering is maximally suppressed. In this case, the DM only scatters off Tl in DAMA. 2 Using the relation σ SI p = 4µ 2 p f 2 p /π, the differential scattering rate can be written as where σ SI p is the proton-DM scattering cross section in the zero momentum transfer limit and µ p is the proton-DM reduced mass.
As discussed above, the inelastic scattering process we consider can be arranged to favor heavier targets. As shown in Eq 2.7 the differential recoil rate is proportional to the proton and nucleon number squared, and the velocity threshold for scattering, v min , decreases with the nucleus mass. For the Tl interpretation of the DAMA signal, the most important effect is the latter. Since Tl is so much heavier than Xe and I, there will be combinations of m χ and δ that give a minimum velocity threshold so large that Xe and I scattering cannot occur, while Tl scattering is still allowed. To illustrate this point, in Fig. 1 we show v min required for 50 keV (left) and 100 keV (right) recoils of W (dotted), Xe (dashed), and Tl (solid), for two different values of m χ , as a function of δ. We see that there exists parameter space in the (m χ , δ)-plane where recoils of a given energy are permitted for the heavier Tl and W targets but kinematically forbidden for the lighter Xe target. Also, the minimum velocity generally increases as the DM mass decreases. Thus, heavier atoms such as Tl will be able to scatter with lighter DM that is kinematically unavailable to lighter atoms like Xe. Note that I and Xe are kinematically similar, and thus the minimum velocity required for DM-I scattering will have the same dependence on the mass splitting as Xe.
To find the differential recoil rate in the detector, we must take into account the finite energy resolution of the experiment. The differential recoil rate in the detector can be rewritten in terms of the electron equivalent of the nuclear recoil energy, E ee , as where φ(E R , E ee ) is the differential response function. Since the DAMA collaboration presents their results corrected for the efficiency of the detector, we have omitted any such correction in our calculation of the recoil rate. The response function is given by where σ E (Q X , E R ) is the energy resolution of the detector and Q X is the quenching factor for a given target nucleus. We take the energy resolution to be with α = (0.448 ± 0.035) √ keV ee and β = (9.1 ± 5.1) × 10 −3 [33]. We use the central values in our calculations and note that changing the resolution has little effect on our conclusions since it is smaller than the size of the energy bins we use after rebinning (see Table 1). However, changing the quenching factor of the detector nuclei has a significant effect on our signal because it causes the recoils to shift from one energy to another.

Quenching Factor
When a DM particle induces a nuclear recoil in a detector, the energy from the recoil is transferred either to electrons or other nuclei. If the energy is transferred to electrons, the recoil can be observed as scintillation light or ionization. If it is transferred to nuclei, the recoil can be observed as phonons or heat. An experiment that measures scintillation light or ionization usually reports the recoil energy in terms of the electron-equivalent energy, E ee , the amount of energy required for an electron recoil to produce the same amount of scintillation light as a nuclear recoil. Generally, nuclear recoils produce less scintillation light than electronic recoils. The true recoil energy and the measured energy are connected by the quenching factor, E ee = Q X E R , which is the most likely value of the measured energy in the distribution corresponding to the response function given by Eq. (2.9). The quenching factor is unique to each nucleus and generally depends on the recoil energy. However, here we assume that the quenching factor is constant as a function of energy for a given type of detector nucleus within the energy range relevant for DAMA.
The Tl quenching factor has never been measured directly. Chang et al. [1] found an approximate range by calculating path lengths and assuming an inverse proportionality to mass. This gives the following relation between the Tl and I quenching factors: 0.62 < Q Tl /Q I < 0.88. Using the range of I quenching factors, 0.06 < Q I < 0.09, we get that 0.037 < Q Tl < 0.079. Given the lack of a direct measurement, we calculate the best-fit The expected Tl signal in the isospin-conserving case for two different quenching factors using the best-fit cross sections. Note that the combination of m χ and δ give a minimum velocity too large for DM-I scattering to occur. The red points with error bars are the DAMA signal in each energy bin. Increasing the quenching factor of Tl corresponds to shifting the zeropoints of the form factor toward larger energies and increasing the period of the oscillatory features in the form factor. This is because a larger Tl quenching factor, Q Tl , for a fixed nuclear recoil energy, E R , corresponds to a larger electron-equivalent energy, E ee = Q Tl E R .
regions for Tl scattering in DAMA with Q Tl ranging from 0.03 to 0.09 with a step of 0.01. As we will demonstrate below, only Q Tl = 0.04 yields a sizeable 2σ best-fit region to the DAMA signal. Thus, we assume this value as a benchmark when investigating the effects of varying other parameters.
In comparison to the recoil spectra of lighter nuclei, Tl recoils are particularly dependent on the quenching factor because of the oscillatory features of the nuclear form factor, F , in the energy range relevant to the signal observed at DAMA. Specifically, the Tl form factor falls to zero in the middle of the DAMA energy range, giving rise to a two-peak shape of the signal spectrum. When the quenching factor increases, the peaks are shifted towards higher energies and become wider. Hence, the Tl quenching factor strongly affects the shape of the expected Tl scattering signal, explaining why the best-fit regions depend so much on the quenching factor. This is illustrated in Fig. 2, where we can see that the same combination of m χ and δ can give two very different signals for two different quenching factors.

Annual Modulation of Dark Matter Velocity Distribution
The velocity distribution of the DM particles in the Earth's rest frame (lab frame), f ( v, t), changes throughout the year due to the orbital motion of the Earth around the Sun. The velocity distribution in the lab frame is found by performing a Galilean boost of the velocity distribution in the DM rest frame,f ( v): is the motion of the lab frame relative to the DM rest frame, v E (t) is the velocity of the Earth relative to the Sun, v = v LSR + v ,pec is the velocity of the Sun relative to the DM rest frame, v LSR = (0, v 0 , 0) is the motion of the Local Standard of Rest (in Galactic coordinates) and v ,pec = (11, 12, 7) km/s is the Sun's peculiar velocity. We take the local circular velocity to be v 0 = 220 km/s. We assume that the DM velocity distribution can be described by a Maxwell-Boltzmann distribution truncated at the DM escape velocity v esc , where the normalization factor is given by Assuming smooth components of the halo, we can divide the differential recoil rate into two terms: One that includes the constant/unmodulated recoil rate, S 0 , and a modulation term, S m , that includes the time-variation of the rate due to the Earth's velocity around the Sun: where t 0 is the date when v E is maximal. The average modulation amplitude in each energy bin of the DAMA detector is then where E 1 and E 2 are the lower and upper energy limits of each bin and the recoil rate is given by Eq. (2.8).
The astrophysical parameters v 0 and v esc , the local circular velocity and the escape velocity of the DM particles, can be subject to significant uncertainties. Typical values used in direct DM detection experiments are v esc = 544 +64 −46 km/s [34,35] and v 0 = 221 ± 18 km/s [36]. As mentioned above we take v 0 = 220 km/s and consider three cases for the maximum DM velocity in the Earth's reference frame denoted as: High, Mid, and Low : v esc + v lab (t 0 ) High: 640km/s + 250km/s Mid: 544km/s + 250km/s Low: 440km/s + 250km/s (2.16) The minimum velocity for DM-Tl scattering given in Eq. 2.1 is dominated by the term δ/ √ 2m X E R when the recoil energy is relatively small, and otherwise depends approximately linearly on E R in the parameter space of interest. This leads to the dependence of v min on the (electron equivalent) recoil energy illustrated in Fig. 3, with a minimum somewhere in the lowenergy part of the spectrum. Scattering is allowed for a given v esc when the minimum velocity is lower than the horizontal lines corresponding to the maximal DM velocity v esc + v lab (t). Lower escape velocities give smaller intervals of E ee where scattering is allowed. Because of the shape of v min , the corresponding recoil spectrum will have less signal in the low (E ee < 2 keV) and high (E ee > 7 keV) energy bins compared to a recoil spectrum with a higher escape velocity. . For scattering to occur, the minimum velocity must be smaller than the maximal DM velocity. When δ increases or m χ decreases, a larger escape velocity is required in order to maximize the signal in each bin. Thus, both the shape and normalization of the signal depend on the minimum velocity. Note that for clarity in these plots we are not accounting for the finite energy resolution of the detector and taking the central value of the energy distribution given by the response function in Eq. (2.9).

Analysis Strategy
In order to improve the signal-to-noise ratio, an alternative binning of the DAMA signal has been used, according to the procedure in Ref. [37]. DAMA presents their recoil spectrum in 36 energy bins, each with a width of 0.5 keV ee . Using 36 bins has two main problems: Most of the bins at high recoil energies have a width smaller than the energy resolution of the detector, and a WIMP signal at DAMA's higher energies is typically negligible. These additional bins thus increase statistical noise without adding information about the signal. The re-binning employed in our analysis groups together adjacent bins with widths substantially more narrow than the energy resolution and groups together all higher-energy bins into one single bin. The resulting bins and their average modulation amplitudes are presented in Table 1.
We have performed a goodness-of-fit test with the theoretical modulation amplitude in Eq. (2.15) and the amplitude detected in the DAMA/LIBRA experiment [38]. We use a χ 2 -test to find the best-fit combinations of m χ , σ SI p , and δ, where S T m,k (m χ , δ, σ SI p ) is the theoretically expected amplitude in bin k, S m,k is the measured amplitude and σ k is the corresponding standard error from the DAMA data. We use the  Table 1: Average modulation amplitudes observed by DAMA over the given energy bins, after rebinning as in Ref. [37]. We have used the modulation amplitudes for the whole data sets: DAMA/NaI, DAMA/LIBRA phase-1 and DAMA/LIBRA phase-2, as presented in Ref. [38]. The rebinning is performed in order to improve the signal-to-noise ratio.
rebinned modulation amplitudes and errors presented in Table 1. For each m χ and δ, we minimize χ 2 with respect to the cross section σ SI p . Note that our confidence regions represent the confidence regions from a perfect fit with the data. The number of degrees of freedom is equal to the number of bins minus the number of free parameters. Using the rebinned data, this gives a total of 7 d.o.f. since we are fitting three parameters: m χ , δ, and σ SI p . We do not include the Tl quenching factor as a free parameter, since we find the confidence regions in (m χ , δ)-space for fixed values of the quenching factor.

Effects of varying astrophysical parameters
In this section, we discuss how varying astrophysical parameters can change the recoil spectra for iDM scattering in the DAMA detector. We focus on the case where isospin-violation suppresses DM-I scattering in order to isolate the effects of the astrophysical parameters on DM-Tl scattering. We keep the local DM density fixed at ρ χ = 0.3 GeV/cm 3 since the associated changes in normalization of recoil spectrum can always be absorbed into the scattering cross section. Note that we only change one parameter at a time, so for all v esc the circular velocity is fixed at the typical value of v 0 = 220 km/s, and when we vary v 0 we keep the escape velocity fixed at v esc = 544 km/s.
As an examination of how well the DAMA signal can be fit in our model, we study the dependence of the best-fit regions on the escape velocity. As shown in the left panel of Fig. 4, the DAMA signal peaks at lower energies and then falls off, with almost no signal and very small errors in the high energy bins. When v esc increases, more signal is allowed in the high energy bins. Since the error is so small in the high energy bins, the fitting procedure will prioritize these data points and choose smaller σ SI p so that the overall normalization of the signal is reduced in order to achieve the best fit. In turn, the smaller normalization yields a worse fit to the low-energy part of the spectrum, where the DAMA signal is largest. As v esc increases, this fit will eventually be so bad that it is no longer within 3σ. The right panel of + +  As one can see from the figure, the fits with higher escape velocities become worse in the low energy bins of the DAMA signal (see the dotted and solid lines corresponding to v esc = 544, 640 km/s), while maintaining a good fit to the high-energy bins. Recalling the earlier discussion of how changes in escape velocity can impact the recoil kinematics for various choices of m χ and δ, we can then analyse the associated changes to the best-fit regions in the (m χ , δ) plane in the right panel of Fig. 4. When the escape velocity increases, the best-fit regions within 3σ increase in size and are shifted towards higher δ. When v esc increases, more regions of parameter space open up for scattering because the threshold, v min < v esc + v lab , allows for higher v min . Since the minimum velocity increases with δ and decreases with m χ , the best-fit regions expand in size towards higher δ and lower DM masses. Due to the increasingly poor fit to the low energy part of the DAMA signal as v esc continues to increase, the best-fit regions also shift towards higher δ and lower m χ when v esc increases.
Similar to v esc , when the circular velocity increases the best-fit regions are expanded to include larger δ and smaller m χ . Since v lab ∼ v E (t) + v 0 , the maximum DM velocity increases when v 0 increases so that more combinations of δ and m χ are allowed. But, since the uncertainties in v 0 are smaller relative to the maximum DM velocity in the Earth's reference frame, this effect is less significant than the corresponding effect for v esc . Also, when v 0 increases the best-fit regions are shifted towards higher δ. This is due to the velocity integral in Eq. (2.7). When v 0 increases the velocity integral, or mean inverse speed, is amplified, but not by the same amount for all energies; the velocity integral is generally more amplified in the high-energy bins than the lower energy bins. This leads to a similar effect as for v esc : when v 0 increases we have too much signal in the high energy bins to find a fit within 3σ to the DAMA data.

Constraints from other direct detection experiments
The lack of a DM signal in all direct detection experiments other than DAMA puts significant constraints on DM models. In this work we have chosen to focus on XENON1T [5,6] and CRESST-II [8] to constrain our best-fit models for DM-Tl inelastic scattering in DAMA. We will show that the lack of a DM signal in these experiments rules out all of the parameter space relevant for DM-Tl inelastic scattering. In this section we present how we have calculated the constraints from XENON1T and CRESST-II, and the effects that varying the astrophysical parameters have on these constraints.
In order to find the regions of (m χ , δ) space where the Tl interpretation of the DAMA signal is ruled out by XENON1T and/or CRESST-II, we calculate the expected number of events for a set of parameters where the Tl interpretation gives a good fit to the DAMA data. The number of expected events is, where dR dEee E ee , t = dR dEee E ee , t, m χ , δ, σ BF p , M T is the exposure of the detector (M is the total target mass and T is the exposure time) and E min and E max are the minimum and maximum recoil energies we consider for a given experiment. For each (m χ , δ) in the parameter space presented in Figs. 6 and 7 we find the best-fit cross section for DM-Tl scattering in DAMA, σ BF p , and calculate the number of expected events in XENON1T and CRESST-II with this cross section. Since the DAMA crystals only contain 0.1% of Tl dopant, the cross sections required to get a good fit to the DAMA data are typically ∼ 1000 times larger than what would be required for DM-I scattering to give a good fit. Using this approach, we find that all the best-fit regions presented in Figs. 6 and 7 are ruled out by a combination of XENON1T and CRESST-II, as they predict a number of events in these experiments that far exceed what has actually been observed.
Below we describe constraints for the isospin conserving and isospin violating benchmarks defined by Eqs. 2.6a and 2.6b. Note we have also considered isospin violating cases tuned such that DM scattering off either Xe or W is maximally suppressed in order relax the respective constraints from XENON1T or CRESST-II. When DM-W scattering is maximally suppressed, the associated regions of best fit to the DAMA signal are similar to those of isospin conserving case described in Sec. 4. When DM-Xe scattering is maximally suppressed, the best fit regions are also similar to those of the isospin conserving case, but the cutoff where DM-I scattering no longer affects the overall signal is at lower m χ and smaller δ. While the associated number of events predicted at XENON1T or CRESST-II in such cases are reduced, the variety of Xe and W isotopes present in either target prevents more than an O(1) suppression. As the number of events predicted at XENON1T and CRESST-II in the best fit regions to the DAMA signal for the isospin conserving case are significantly larger llll Figure 5: The regions in the (m χ , δ)-plane excluded by XENON1T [6,39] (left) and the remaining relevant regions excluded by CRESST-II [8] (right) for different v esc , assuming isospin-conserving DM. The maximum velocity of the DM particles is linearly dependent on v esc . When we increase v esc , regions of parameter space previously forbidden by the threshold v min < v esc + v lab (t) for DM-nucleus scattering become available. For the case of isospinviolating DM, the regions constrained by XENON1T and CRESST-II are approximately the same as what is presented here.
than what are observed, all such points remain excluded even when DM scattering off Xe or W is maximally suppressed.

XENON1T
In order to constrain nuclear recoils over a larger energy range than considered in Ref. [6], we use data from Ref. [39] with an accumulated exposure of 0.65 ton×years. Although the latter analysis shows an excess of electron recoil events, virtually no events arising from nuclear recoils are detected. To set conservative bounds we assume that approximately 1% of the ∼ 42000 electron recoils in Ref. [39] could be misidentified as nuclear recoils [40]. Given a best-fit cross section for DM-Tl scattering at DAMA, the associated combination of (m χ , δ) is ruled out by XENON1T if the number of expected events corresponds to more than a 3σ statistical fluctuation in the number of misidentified electron recoils, approximately 500 events.
Although we will show in Section 4 that the combined constraints from XENON1T and CRESST-II rule out all of the best-fit regions to the DAMA signal, here we focus on the points within the best-fit regions of (m χ , δ)-space which are robustly constrained by XENON1T alone. When considering the large cross sections which yield the best fit to the DAMA signal in combination with the large exposure of XENON1T, such points ruled out by XENON1T alone typically predict at least 1000 events. For these points, we have checked that allowing for the cross sections to be reduced by 20% while still keeping a 3σ fit to the DAMA signal results in a prediction of at least 500 events at XENON1T and are thus still ruled out before considering the constraints from CRESST-II. Near the regions of (m χ , δ)-space where DM-Xe scattering is kinematically forbidden, the expected signal in XENON1T falls off rapidly. Thus, the parts of parameter space that are not ruled out by XENON1T are close to where the kinematics of inelastic DM-Xe scattering do not allow for recoil energies of 1 − 201 keV.
As discussed in Section 2.2, the DM escape velocity and the local circular velocity are both subject to significant uncertainties (see Refs. [41,42] for recent analyses of astrophysical uncertainties in XENON1T constraints). The regions of parameter space that are ruled out by other direct detection experiments also depend on these parameters. In particular, the constraints from XENON1T depend heavily on the choice of astrophysical parameters since the excluded regions are determined almost entirely by kinematics. As shown in the left panel of Fig. 5, regions of the (m χ ,δ) plane which are excluded by XENON1T grow when the escape velocity increases (there are similar effects when increasing the local circular velocity). This is a kinematical effect: The XENON1T constraints depend on the threshold for scattering, given by v esc +v lab (t) > v min . According to Eq. 2.1, the minimum velocity increases with δ and decreases with m χ . Thus, when v 0 and v esc increase, the regions constrained by XENON1T increase in size and include larger δ and smaller m χ .

CRESST-II
We follow a similar procedure to calculate the constraints from CRESST-II. 3 We use the data from the TUM-40 detector in CRESST-II because it has the fewest observed events in the acceptance region for dark matter scattering (E R = 0.6 − 40keV) and the highest exposure after cuts have been applied. We use the data from Ref. [8,43] with an exposure of 29.35 kg × days, and a cut-survival probability of 0.8 for DM-Tungsten (W) scattering. There is an excess of events in CRESST-II for recoil energies E R 1 keV, but these cannot be distinguished from the e − /γ-background. To make sure no events within the 90% confidence region of the e − /γ-band are included, we only consider nuclear recoils with E R 10 keV. There are no events above this limit within the signal acceptance region of the CRESST-II analysis. Given the associated best-fit cross section for DM-Tl scattering in DAMA at a point in (m χ , δ) space, we conservatively demand that the expected number of events in CRESST-II is less than 10 events.
We find that CRESST-II rules out the remaining (m χ , δ)-points in the best-fit regions, discussed in detail in Sec. 4, which are not excluded by XENON1T alone. The smallest number of events predicted in CRESST-II at such points is 40 assuming the best-fit cross sections to the DAMA signal at each point. The cross sections can be reduced by 20% while still keeping a 3σ fit to the DAMA signal at each point and the smallest number of predicted events decreases to 34. Thus, the combined constraints from XENON1T and CRESST-II on the best-fit regions of (m χ , δ)-space are robust to any potential reduction of the cross sections from the best-fit values while still maintaining a 3σ fit to the DAMA signal.
The constraints from CRESST-II are shown in the right panel of Fig. 5. The regions constrained by CRESST-II behave in a similar manner to those from XENON1T when v esc and/or v 0 increases. In general, the constraints from CRESST-II are less robust than those from XENON1T since CRESST-II has a much smaller exposure. However, the CRESST-II constraints are particularly important in the region of parameter space where XENON1T does not rule out the best-fit regions entirely. More specifically, when Q Tl = 0.05 − 0.09 XENON1T cannot entirely rule out the the associated best-fit regions due to kinematic effects. However, CRESST-II contains W target nuclei, which are more similar in mass to Tl than Xe. Thus, any best-fit regions not excluded by XENON1T are excluded by CRESST-II since W has a large enough mass for scattering with DM to occur. We also note that the regions of (m χ , δ)-space not ruled out by CRESST-II are somewhat less sensitive to the kinematics of inelastic DM-W scattering and are more a function of the number of events decreasing with the best-fit cross section to the DAMA signal at lower m χ and higher δ.

Best-fit regions for DAMA modulation signal
The results presented in the current section are based on the expected annual modulation signal in DAMA for v esc = 544 km/s (Mid) and 640 km/s (High), v 0 = 220 km/s and ρ χ = 0.3 GeV/cm 3 . We identify the best-fit regions for two cases: isospin-conserving DM with f n /f p = 1, and isospin-violating DM with f n /f p = −53 127−53 . For the isospin violating case, f n /f p is chosen so that the effective coupling to I vanishes and, thus, the only nonnegligible contribution in DAMA comes from DM-Tl scattering. While we allow for DM scattering off of I in the isospin-conserving case, we focus on the parameter space of the iDM model where the dominant scatting contribution comes from DM-Tl scattering. Because the scattering kinematics are very similar for I and Xe target nuclei, fits to the DAMA signal arising primarily from scattering off I are rather trivially ruled out by the XENON1T constraints described in the previous section. The best-fit points for each value of Q Tl are given in Table 2.
In Figure 6 we show the only interesting 2σ best-fit region within the various scenarios we have explored, the isospin-violating case with Q Tl = 0.04. For any of the other Tl quenching factors we consider here, there are no substantial regions of parameter space in which the dominant scattering contribution comes from DM-Tl scattering that fit the DAMA annual modulation signal at 2σ. 4 When Q Tl = 0.05, there is a tiny 2σ region in the isospinconserving case. This region is a thin line centered near m χ 80 GeV and δ 171 keV, where the parameter space is ruled out by CRESST-II. 5 For both escape velocities considered in our analysis, we see that the 2σ best-fit regions in the isospin-violating case with Q Tl = 0.04 have 165 keV δ 190 keV (up to 195 keV when v esc = 640 km/s) and extend indefinitely towards higher m χ . As m χ increases, the expected DM-Tl scattering signal reaches a plateau and no longer changes with increasing DM mass. However, the 2σ best-fit regions for Q Tl = 0.04 and both escape velocities are completely ruled out by XENON1T, as the suppression of scattering off of Xe targets becomes less effective for larger DM masses.
The best-fit regions within 3σ are presented in Fig. 7. The 3σ best-fit regions for the isospin-violating case with Q Tl = 0.04 extend to significantly lower DM masses, m χ 100 GeV, than the corresponding 2σ best-fit regions. Also, at even lower DM masses, there are several smaller 3σ best-fit regions with Q Tl > 0.04. In isospin-conserving cases, we see that the 3σ best-fit regions vanish for smaller quenching factors because DM-I scattering becomes kinematically allowed at higher DM masses but does not yield a good fit to the DAMA modulation signal.  Comparing Figs. 5 and 7, we see that XENON1T rules out all of the 3σ best-fit regions for m χ 100 GeV and part of the 3σ best-fit regions for m χ 100 GeV, while CRESST-II rules out all lower mass 3σ best-fit regions. Since CRESST-II contains W, which has mass nearly as large as Tl, DM-W scattering is kinematically allowed for approximately all mass and δ ranges where DM-Tl scattering occurs. Thus, the parts of the 3σ best-fit regions not excluded by XENON1T are excluded by CRESST-II. We therefore find that a combination of XENON and CRESST data rules out all regions of parameter space for which inelastic scattering could explain the DAMA data.

Conclusions
In this paper we have studied the compatibility of the observed DAMA modulation signal with inelastic dark matter (DM) scattering in DAMA. This type of scattering has been proposed Figure 6: The best-fit regions within 2σ for DM-Tl scattering where DM-I scattering is maximally suppressed by isospin violation, as described in Eq. (2.6b). Left: The bestfit region assuming the astrophysical parameters ρ χ = 0.3 GeV/cm 3 , v 0 = 220 km/s and v esc = 544 km/s. Right: The best-fit region assuming v esc = 640 km/s and keeping all other astrophysical parameters the same. In both cases, there are no substantial regions of parameter space within 2σ for Tl quenching factors other than Q Tl = 0.04. The best-fit regions extend indefinitely towards larger m χ , with the corresponding range of δ remaining constant. The shaded regions in both panels, corresponding to the only substantial iDM parameter space where DM-Tl scattering fits the DAMA signal within 2σ, are ruled out by the constraints from XENON1T shown in Fig. 5.
as a solution to the discrepancy between the observed signal in DAMA and the lack of signal in other direct detection experiments. In particular, it has been proposed that the DM primarily scatters off of the 0.1% Thallium (Tl) dopant in the DAMA crystals. Since inelastic scattering prefers heavy nuclei, it is possible for DM scattering off of atoms much lighter than Tl, such as xenon, to be kinematically forbidden. Thus, if inelastic scattering gives the leading contribution to the DM-nucleon scattering cross-section, the constraints from liquid noble gas target experiments, such as XENON1T [5,6], can be suppressed. Previously, Chang, Lang and Weiner [1] in 2010 and more recently the DAMA collaboration [2] in 2019 found regions in parameter space where DAMA data was consistent with an interpretation of inelastic scattering of DM with Tl in the detector. In light of the recent DAMA results, we have expanded upon these previous analyses by using the most recent data from DAMA/LIBRA phase I and II [12,44] and including constraints from both XENON1T [6,39] and CRESST-II [8].
We have compared the expected modulation signal for inelastic DM scattering in DAMA with the observed signal for two cases: Isospin-conserving DM and isospin-violating DM, where DM-I scattering is maximally suppressed. The main sources of uncertainty are the quenching factor of Tl and the astrophysical parameters v 0 , v esc , and ρ χ . The quenching factor of Tl has never been measured directly, but an approximate range was found by Chang et al. Figure 7: The best-fit regions within 3σ for isospin violating (top) and isospin conserving (bottom) scattering in DAMA, assuming either v esc = 544 km/s (left) or v esc = 640 km/s (right). Note that the 3σ regions for Q Tl = 0.04 extend indefinitely towards higher m χ in the same manner as the 2σ regions presented in Fig. 6, with the same range of δ at larger m χ . Note there are no best-fit regions within 3σ when Q Tl = 0.02 and 0.03. For the isospin conserving case, the regions within 3σ corresponding to different Tl quenching factors come from DM-Tl scattering only since the minimum velocities required to scatter in these parts of parameter space are too large for DM-I scattering to occur. This is illustrated by the sharp cutoff starting at m χ = 70 GeV and δ ∼ 160 GeV in the lower left panel (starting at m χ ∼ 60 GeV and δ ∼ 175 keV in the lower right panel) which shows where the DM-I scattering becomes kinematically forbidden. All regions shown in this figure are ruled out by a combination of the XENON1T and CRESST-II constraints shown in Fig. 5. in their 2010 paper. As a result of their calculations, we have performed our analysis for a range of Tl quenching factors from 0.03 to 0.09 with a step size of 0.01. We have shown that the best-fit regions depend heavily on the choice of quenching factor, since the quenching factor directly affects the shape of the modulation signal (see Fig. 2). We have also considered how the fits to the DAMA signal depend on the astrophysical parameters v esc and v 0 . We have found that the best-fit regions increase slightly in size and are shifted towards higher mass splitting and lower DM masses when these two parameters are increased.
We find that all regions in (m χ , δ) space that give a good fit to the DAMA modulation signal are ruled out by XENON1T, CRESST-II, or both. For the isospin-violating case there is a 2σ region when Q Tl = 0.04. For the isospin-conserving case, we have found a 2σ region for δ < 120 keV that comes from DM-I scattering. Constraints from XENON1T exclude all of the 2σ best-fit regions. The 3σ best-fit regions are ruled out by XENON1T (high mass) and/or CRESST-II (low mass). The constraints from XENON1T and CRESST-II have been found by directly comparing the number of events that have been observed in these experiments with what is expected given the Tl interpretation of the DAMA signal. The best-fit regions presented in the paper by the DAMA collaboration [2] are also within the parameter space that is ruled out by XENON1T and/or CRESST-II. Based on these results, we conclude that inelastic DM-Tl scattering in DAMA cannot explain the discrepancy between the modulation signal in DAMA and the null results of other experiments.