Strongly interacting dark matter and the DAMA signal

TL;DR

This paper demonstrates that strongly interacting dark matter models designed to explain the DAMA signal cannot reproduce its observed time dependence, due to constraints from detection experiments and diffusion effects.

Contribution

It provides a detailed analysis showing the incompatibility of SIMP dark matter models with the DAMA signal's temporal modulation features.

Findings

Models predict incorrect phase of annual modulation.

Models predict larger diurnal modulation amplitude.

Constraints from direct detection limit model parameters.

Abstract

We show that models of strongly interacting (SIMP) dark matter built to reproduce the DAMA signal actually cannot account for its time dependence. We discuss the constraints on this type of models coming from direct detection experiments and study the propagation of thermalised dark matter particles in the ground for the allowed values of the parameters. We consider a simple 1D diffusion and a more detailed 3D diffusion. In both cases the predicted signal has either the wrong phase of the annual modulation or a much larger amplitude of the diurnal modulation.

Equations

$$l_{\mathrm{th}}=\left(\frac{M}{m}\right)\log \left(\frac{v_0}{v_{\mathrm{th}}}\right)\frac{1}{n{\sigma }_N},$$

$$\frac{\partial N(\vec{x},t)}{\partial t}=D\Delta N(\vec{x},t)-v_{\mathrm{d}}\frac{\partial N(\vec{x},t)}{\partial z}+f(\vec{x},t),$$

$$D=\frac{\pi }{8}\frac{M}{\mu }\lambda v_{\mathrm{th}}=\frac{m+M}{mn{\sigma }_{\mathrm{N}}}\sqrt{\frac{\pi kT}{8M}},$$

$$v_{\mathrm{d}}=\frac{Mg\lambda }{\mu v_{\mathrm{th}}}.$$

$$f(\vec{x},t)=n_{\mathrm{loc}}\iint d{x}^{\ast }d{y}^{\ast }{\omega }_l(\overrightarrow{x^{\ast }}-\vec{x},t)v\left(\left|\overrightarrow{x^{\ast }}-\vec{x}\right|\right).$$

$$N(z,t)={\int }_{-\infty }^{t}d\tau {\int }_0^{\infty }d\xi \exp \left(\frac{2v_{\mathrm{d}}(z-\xi )(t-\tau )-{v_{\mathrm{d}}}^2(t-\tau {)}^2}{4D(t-\tau )}\right)G(z-\xi ,t-\tau )f(\xi ,\tau ).$$

$$G(z-\xi ,t-\tau )=\frac{1}{\sqrt{4\pi D(t-\tau )}}\left[\exp \left(-\frac{(z-\xi {)}^2}{4D(t-\tau )}\right)-\exp \left(-\frac{(z+\xi {)}^2}{4D(t-\tau )}\right)\right],$$

$$G(x-\chi ,y-\upsilon ,z-\xi ,t-\tau )=\frac{1}{(4\pi D(t-\tau ){)}^{3/2}}\exp \left(\frac{(x-\chi {)}^2+(y-\upsilon {)}^2+(z-\xi {)}^2}{4D(t-\tau )}\right).$$

$$f(\chi ,\upsilon ,\xi ,t)\approx n_{\mathrm{loc}}S\sum _i\sum _j{\omega }_l(x_i-\chi ,y_j-\upsilon ,z_{ij}-\xi ,t)v(l),$$

$$\begin{gathered} N(x,y,z,t)\approx n_{\mathrm{loc}}S\sum _i\sum _j{\int }_{-\infty }^{t}d\tau {\int }_{-\infty }^{\infty }d\chi {\int }_{-\infty }^{\infty }d\upsilon {\int }_0^{\infty }d\xi \times \\ \times G(x-\chi ,y-\upsilon ,z-\xi ,t-\tau )\exp \left(\frac{2v_{\mathrm{d}}(z-\xi )(t-\tau )-{v_{\mathrm{d}}}^2(t-\tau {)}^2}{4D(t-\tau )}\right)\times \\ \times {\omega }_l(x_i-\chi ,y_j-\upsilon ,z_{ij}-\xi ,t)v(l). \end{gathered}$$

$$\frac{\partial v}{\partial t}=D\frac{{\partial }^{2}v}{\partial z^2}+\beta \frac{\partial v}{\partial z}.$$

$$\frac{\partial u}{\partial t}=D\frac{{\partial }^{2}u}{\partial z^2},$$

$$v=\exp (\mu z+\lambda t)u\text{ with }\mu =-\frac{\beta }{2D}\text{ and }\lambda =-\frac{{\beta }^2}{4D}.$$

$$L^{\prime }=\frac{\partial }{\partial t}-D\frac{{\partial }^{2}v}{\partial z^2}-\beta \frac{\partial v}{\partial z},$$

$$L=\frac{\partial }{\partial t}-D\frac{{\partial }^2}{\partial z^2}.$$

$$G^{\prime }=\exp (\mu (z-\xi )+\lambda (t-\tau ))G.$$

$$L^{\prime }{G}^{\prime }=\exp (\mu (z-\xi )+\lambda (t-\tau ))\delta (z-\xi )\delta (t-\tau ),$$

$$\frac{\partial v}{\partial t}=D\frac{{\partial }^{2}v}{\partial z^2}+\beta \frac{\partial v}{\partial z}+f,$$

$$v(z,t)=\int d\tau \int d\xi \exp \left(-\frac{\beta (z-\xi )}{2D}-\frac{{\beta }^2(t-\tau )}{4D}\right)G(z,\xi ,t,\tau )f(\xi ,\tau ),$$

$$J^{+}=\frac{⟨N^{+}/A⟩}{⟨\tau ⟩}.$$

$$\frac{1}{2}d\cos \theta dln(l),$$

$$P(S)=\frac{1}{\lambda }\exp \left(-\frac{S}{\lambda }\right),$$

$${\int }_{l/\cos \theta }^{\infty }dS\frac{1}{\lambda }\exp \left(-\frac{S}{\lambda }\right)=\exp \left(-\frac{l}{\lambda \cos \theta }\right).$$

$$\lambda +q\lambda +q^2\lambda +\dots =\lambda (1+q+q^2+\dots )=\frac{\lambda }{1-q}=k\lambda ,$$

$$⟨N^{+}/A⟩=\frac{1}{2}{\int }_0^{1}d\cos \theta {\int }_0^{\infty }dl\exp \left(-\frac{l}{k\lambda \cos \theta }\right)n(l).$$

$$n(l)\approx n(z_0)-\frac{\partial n}{\partial z}l,$$

$$⟨N^{+}/A⟩=\frac{k\lambda }{2}\left(\frac{n(z_0)}{2}-\frac{k\lambda }{3}\frac{\partial n}{\partial z}\right).$$

$$⟨\tau ⟩=\int dSP(S)\int dvP(v)\frac{S}{v}=\frac{4k\lambda }{\pi ⟨v⟩},$$

$$J^{+}=\frac{\pi ⟨v⟩}{8}\left(\frac{n(z_0)}{2}-\frac{k\lambda }{3}\frac{\partial n}{\partial z}\right).$$

$$J=J^{+}-J^{-}=-\frac{\pi k\lambda ⟨v⟩}{12}\frac{\partial n}{\partial z}.$$

Full text

1 Introduction

The DAMA and DAMA/LIBRA results have been with us for 20 years [1], and at present the observed modulation is the most significant unexplained experimental result in physics, reaching a level of 12.9 $\sigma$ [2]. The signal consists of hits in NaI scintillators, and corresponds to the emission of electromagnetic radiation with energies in the range 1 to 6 keV. To disentangle this signal from backgrounds, DAMA/LIBRA monitors it as a function of time. If the signal is due to weakly interacting massive particles (WIMP) scattering on the detector nuclei, and if dark matter is at rest with respect to our galaxy, then the motion of Earth in the galaxy creates a specific time signature, which is the signature. More specifically, the Sun moves towards the Cygnus constellation at roughly 220 km/s, and the Earth moves around the Sun at about 30 km/s. Depending on the time of the year, these velocities will add or subtract. As the Earth velocity is never exactly parallel or antiparallel to that of the Sun, only the projection of the velocity onto the Sun’s velocity modulates the flux, and this amounts to a variation of about 6 % of the flux of dark matter. What DAMA has been detecting for 20 years is such a periodic signal, and it has precisely the correct phase, being maximum around June 2nd and minimum around December 2nd. Besides this modulation, the bulk of the dark-matter flux and the residual backgrounds should also produce events.

This could have been the answer to the 100-year old puzzle of dark matter, but the rub has been that no other detector sees a signal, and several of them are severely incompatible with the DAMA/LIBRA results111Some studies indicate that the new DAMA data [2] is itself inconsistent with the conventional isospin-conserving spin-independent WIMP interaction with nuclei [3]. [4, 5, 6]. A number of explanations of this problem have been tried, assuming that the DAMA signal is due to conventional physics, but all of them have been rebutted [2, 7]. As no standard physics explanation seemed to work, several groups proposed models that could lead to a signal in DAMA and not in other detectors. We are aware of three classes of such models: mirror matter [8, 9], resonant dark matter [10], and specific strongly interacting massive particle (SIMP) models [11, 12]. In this paper, we shall concentrate on the latter.

Models based on SIMPs find their origin in the work of M. Khlopov and collaborators, who proposed that a composite bound state, made of a SIMP of charge $-2$ and a 4He nucleus, could be the dominant form of dark matter [13] and the key to the explanation of the DAMA signal [14, 15]. The elastic cross section of this OHe atom is comparable to that of neutrons, and after hitting the ground it will quickly thermalise and turn into a cloud of slowly-moving heavy particles. As it is moving with thermal energy, it does not create detectable recoils when it hits nuclei. Instead, it was argued [14] that a repulsive force would arise when OHe got close to a nucleus, and that it might allow bound states with sodium. The formation of these bound states would result in the emission of photons of a few keV energy, which would make the DAMA signal. The demise of this model came from the absence of such a repulsive force [16], which would have observable and disastrous consequences [17].

But the essential features of this model could be preserved, in the context of a dark sector, with a subdominant component made of dark atoms [11] or dark antiatoms [12]. The elastic cross sections could be made large enough for the incoming composite SIMPs to thermalise and escape detection from nuclear recoil. It was also possible to produce a signal via binding with some of the elements of the DAMA detector, and not with those of other detectors. These two models are examples of a class of models which could, in principle, reproduce the DAMA signal. However, as we shall show in this paper, the first ingredient, i.e. the thermalisation of the SIMPs before they reach the detector, makes it impossible to reproduce the time signature observed by DAMA [18].

This paper is organised as follows: in Section 2, we spell out the constraints on the elastic cross section of SIMPs, in Section 3 we explain the various propagation models that we used, and show the time dependence of the signal according to the degree of precision of the diffusion model, in Section 4 we discuss the impact of the gravitational focusing of the Sun on the signal. We then conclude about the consequences of this work for SIMPs and composite dark matter.

2 Constraints

As the problem will come from thermalisation, we can study a generic class of models. We do not need to discuss the particular mechanism which creates the signal in the detector and we can consider only the parameters that determine the propagation, i.e. the particle mass $M$ and the cross section $\sigma_{N}$ of elastic scattering on the nuclei of the ground. The depth of thermalisation is given by

[TABLE]

with $m$ the mass of atoms in the crust, $v_{0}$ the incoming velocity, $v_{th}$ the thermal velocity, and $n$ the density of the crust. For example, to thermalise at the depth of LNGS (1.4 km) a particle with a mass of 1 TeV should have $\sigma\approx 5\ 10^{-26}{\,\rm cm}^{2}$. The first set of constraints, shown in Fig. 1, actually comes from underground direct-detection experiments with null results. The experiments that are the most sensitive222This is due to a combination of multiple factors, such as the threshold energy, detector material and the depth of the laboratory. to the particles we consider are XENON100 ($E_{\rm th}=700{\,\rm eV}$) [21] and CRESST-III ($E_{\rm th}\approx 100{\,\rm eV}$) [20], located at the same site as DAMA.

Experiments operating at the top of the Earth atmosphere or above provide the second set of constraint on SIMPs, and limit the parameter space from above. The most severe constraint comes from the balloon-borne direct detection experiment RRS [19]. For example, the aforementioned cross section of $\sim 10^{-26}{\,\rm cm}^{2}$ lies a few orders of magnitude above the limit for 1 TeV particle. To loosen this constraint one has to assume that SIMPs are a subdominant fraction of the local dark-matter density. Taking into account the first set of constraints this fraction has to be $<5\%$. Hereafter we take the fraction of SIMPs to be 1%, although we shall also study one particular case where this fraction is $0.1\%$. Note, that these constraints are for the elastic scattering cross section on silicon nuclei and were taken from the analysis of the RRS data [19], not recalculated from the constraints on dark-matter-nucleon cross section (e.g., as in [22]). Since it is a reasonable assumption to take the average rock as composed entirely of silicon atoms (see below), we can apply these constraints directly.

3 Propagation models

3.1 Ray propagation

To explain the problem, let us first use a very simplified picture of propagation. Assume that the ground can be treated as an infinite plane, on which an infinite beam of dark matter falls, and assume that all the particles have the same velocity $V_{\rm DM}=v_{0}$. This beam then hits the ground, the particles slow down, and after a while move down because of gravity. Consider a packet of SIMPs travelling through the ground, as in Fig. 2. As the speed is lower below ground, the packet gets concentrated. But the packet also gets spread, as its section is inversely proportional to $\cos\theta$, with $\theta$ the angle with the azimuth, and becomes zero if $|\theta|>\pi/2$, i.e. when the flux of dark matter comes from below the horizon. So the density of the packet is proportional to $\cos\theta$. $\theta$ is a complicated function of time which can be found in [18]. This function depends on the Earth rotation, and the flux of dark matter is screened by the Earth part of the day for a period of the year. This leads to a daily variation, and to a phase shift in the yearly flux.

The resulting time profile of the expected signal in the DAMA detector is shown in Fig. 3, where we compare it to the WIMP signal and to the results of the diffusion approximation, which we shall discuss later. One immediately notices a substantial time delay of the annual modulation (left plot) compared to the WIMP case that reproduces the DAMA data [23], as well as a non negligible diurnal modulation of the signal (right plot). A modulation of the daily signal is actually present in the WIMP case as well, because of the change of flux due to the Earth rotation, but its amplitude is tiny.

We treat this result as an example of the difficulties that a more detailed model of propagation has to overcome in order to provide a better fit of the DAMA data. One could hope that diffusion will spread the daily variation enough to make it unobservable. The annual phase shift could also be affected by a more precise treatment, which we shall consider now.

3.2 Diffusion

A more realistic treatment of the propagation of thermalised SIMPs in the ground has to account for the fact that it takes time for these particles to reach the detector from the point where they are in equilibrium with the surrounding matter. Furthermore, they constantly scatter on nuclei in the medium and get driven by the gravitational field, so they do not proceed in the same direction as the incoming flux of dark matter. In principle, a fraction of the incoming flux can even scatter back up, out of the ground. We shall assume that once a SIMP gets above the ground it is never coming back (at least, in the same vicinity), so particles, which thermalise close to the surface have lower chances to contribute to the signal.

The processes, described above, resemble very much the diffusion of a gas in the presence of a gravitational field, so we are going to adopt the corresponding diffusion equation for the number density of SIMPs as a function of space and time $N(\vec{x},t)$ to see how it evolves inside the detector:

[TABLE]

where the first term on the right-hand side accounts for diffusion, with $D$ the diffusion coefficient, the second term describes the influence of gravity, with ${v_{\rm d}}$ the drift velocity (parallel to the z-axis) and the third term $f(\vec{x},t)$ is the source function. One should, in principle, solve Eq. (2) in an inhomogeneous unbounded space, as the air and the ground have different diffusion coefficients, for a source function steady in time333The solution of the heat transfer equation in a medium, consisting of regions with different values of thermal conductivity, was found by Sommerfeld in [24].. However, the density of air is much smaller than that of rock, so we can neglect the propagation in the atmosphere and solve the diffusion equation in the space limited by the surface $S_{\rm b}(x,y)$ between the two media with the condition $\left.N(\vec{x},t)\right|_{S_{\rm b}}=0$ imposed on this boundary. Physically, this conditions corresponds to a sink of particles at the boundary.

The diffusion coefficient depends on several parameters: those characterizing the medium, the mass of the SIMPs and the cross section describing their interaction with nuclei in the medium. We show in Appendix B that the following formula holds:

[TABLE]

where $\mu=mM/(m+M)$ is the reduced mass, $n$ is the number density of the medium, $\lambda=(n\sigma_{\rm N})^{-1}$ is the mean free path and $v_{\rm th}$ is the average thermal velocity of SIMPs. Using the Einstein relation [25] $D=\alpha kT$, where $\alpha$ is the ratio of the drift velocity ${v_{\rm d}}$ to the applied force $F=Mg$ with $g=9.81$ m/s2 the gravitational acceleration, one gets the expression of the drift velocity in the gravitational field

[TABLE]

For simplicity, we assume that the rock is composed only of silicon atoms with $m=28{\,\rm GeV}$, following [26]. Although oxygen atoms are actually about 4 times more abundant in rock than silicon atoms, we have checked that our results are insensitive to this, and that our main conclusions remain the same.

The source function $f(\vec{x},t)$ describes the rate of the thermalised SIMPs density increment and depends on many factors. First of all, it depends on the local number density of SIMPs, which we take to be $n_{\rm loc}=0.01\cdot\rho_{\rm loc}/M$, as explained in Sec. 2, where $\rho_{\rm loc}=0.39{\,\rm GeV}/{\,\rm cm}^{3}$ is the conventional value of the local dark-matter density. Second, SIMPs are coming with different velocities, which follow the velocity distribution $\omega_{v}$, hence the distance they travel through the ground until thermalisation also differs. Let $\vec{l}$ be the vector connecting the point $(x^{*},y^{*})$ on the surface $S_{\rm b}$, where the particle enters the ground, with the point $(x,y)$, where it acquires the thermal velocity. It is convenient to transform the velocity distribution into the path length distribution $\omega_{v}\rightarrow\omega_{l}$ (see Appendix C). Third, the rate of SIMPs crossing the surface $S_{\rm b}$ at the point $\left\{x^{*},y^{*}\right\}$ depends on the angle between the velocity of the incoming beam of SIMPs and the normal vector to the surface at the given point $\vec{n}$.

The source function can be constructed in the following way

[TABLE]

The complexity of the solution of Eq. (2) depends mainly on the shape of the boundary surface $S_{b}$. We are going to consider two relevant cases.

3.2.1 The plateau approximation

The simplest shape of the boundary surface $S_{b}$ is obviously a plane, for example $z=0$. Since the total incident flux of SIMPs through the boundary is uniformly distributed on that surface the diffusion is symmetric in $x$ and $y$ and the density gradient is directed along the $z$-axis. In this case Eq. (2) becomes one-dimensional and its solution has the following form

[TABLE]

Here $G(z-\xi,t-\tau)$ is the Green’s function of the one-dimensional diffusion equation in the semi-infinite space with the boundary condition $N(0,t)=0$ [27, p. 209]

[TABLE]

which describes the propagation of particles created at the point $\xi$ at the moment $\tau$. The income of thermalised SIMPs is determined by the source function (5), which also becomes effectively one-dimensional. The factor in front of the Green’s function in Eq. (6) appears due to the presence of the drift term in the Eq. (2) (see Appendix A).

The integration in Eq. (6) can be performed numerically. As an example, we demonstrate the time dependence of the particle density in the DAMA detector ($z=1.4{\,\rm km}$) for SIMPs with $M=1{\,\rm TeV}$ and $\sigma=60{\,\rm mb}$ (Fig. 4). We assume that the signal in DAMA is proportional to the density of SIMPs in the detector $S=\beta N$, where $\beta$ is an additional parameter, which depends on the underlying physics of SIMP-nucleus inelastic interactions. In our analysis we only use this parameter to fit the predicted signal to the annual modulation data and we do not impose any constraints on its value. Following the procedure described in [18] we calculate the residual rate of events and compare it to the DAMA data on annual and diurnal modulation in $2-6$ keV energy interval.

Fig. 5 shows that for $M=1{\,\rm TeV}$ and $\sigma=60{\,\rm mb}$ not only the amplitude of the diurnal modulation of the signal exceeds the data points by orders of magnitude, but also the phase of the annual modulation of the signal is off the best-fit value by about 100 days. The daily variation comes from the fact that the density of dark matter hitting the ground depends on the angle between the flux of dark matter and the ground. The phase of the yearly signal comes from a change in the daily average of that angle due to the obliquity of the ecliptic. Diffusion of dark matter is not strong enough to flatten the first variation, and it cannot change the second. This is mainly because the value of the diffusion coefficient allowed by the constraints in Fig. 1 is quite large $D\gtrsim 1000{\,\rm m}^{2}/{\,\rm s}$, which corresponds to ${v_{\rm d}}\sim 10{\,\rm m}/{\,\rm s}$. Thus, it takes a cloud of SIMPs around a few minutes to get from the surface to the laboratory, which is totally negligible for the relevant observational time scales.

Qualitatively similar results are obtained for the whole allowed range of parameters above the black dashed line in Fig. 1. Returning to Fig. 3, we see that the profile of the signal, predicted within the diffusion approach, does not substantially differ from the much simpler model that we considered before.

As the cross section gets smaller and the thermalisation depth distribution extends much deeper than DAMA, the diurnal modulation of the signal almost vanishes, though the phase of the annual modulation is still wrong (see Fig. 6). Furthermore, the average density of SIMPs in the detector in this case is many orders of magnitude smaller than for $\sigma=60{\,\rm mb}$. Hence, to account for the observed rate, some extreme values of the parameter $\beta$ are required, which, in its turn, implies unrealistically high values of the inelastic cross section. Reducing the fraction of SIMPs to $0.1\%$ to lift the RRS constraint and allow larger cross sections neither recovers the right phase of the annual modulation, nor eliminates the diurnal modulation.

3.2.2 Realistic Gran Sasso terrain surface

Clearly, the surface of the Earth around LNGS is far from being a simple plane. Finding the exact analytical solution (especially, the Green’s function) of the diffusion equation (2) for a realistic humpy terrain as the one surrounding LNGS (see Fig. 7) is very complicated. The problem can be drastically simplified if one neglects the loss of particles through the boundary into the atmosphere444We have studied this simplification in the plateau geometry and our results indicate that the difference between the exact solution and the simplified one is not larger than $\approx 10\%$. and uses a simple Green’s function of the 3D diffusion in the unbounded space

[TABLE]

To further simplify the calculation we break the considered landscape into squares of the size $S=500\times 500{\rm\,m}^{2}$ and construct the source function piecewise

[TABLE]

where $x_{i}$ and $y_{j}$ are the coordinates of the centre of each square and $z_{ij}$ is their elevation. The density of SIMPs at the given point in time and space can be calculated as follows555For the sake of brevity, in this expression we do not explicitly indicate that the source distribution is truncated above the boundary surface.

[TABLE]

Following the same routine as described above we compare our predictions to the DAMA data. The results we get for the realistic surface are not very different from those obtained in the plateau approximation (see Figs. 8 and 9). Thus, a more detailed approach doesn’t solve the aforementioned problems with the DAMA explanation.

4 Gravitational focusing

Finally, we have considered the possibility that the annual phase is modified because the dark-matter flux must be bent by the Earth and the Sun (see e.g. [29, 30]). The magnitude of this effect is inversely proportional to the square of the velocity that an incoming particle has far away from the Sun and the maximum is expected around March. Since the maximum of the unperturbed annual dark matter velocity modulation is at the beginning of June, the effect of GF can modify the phase of the event rate in the direct-detection experiment like DAMA depending on the mass of dark-matter particles. In a WIMP scenario, this sets the bound on the velocity that a WIMP can have in order to create a recoil at the threshold energy. Unlike WIMPs, the signal from SIMPs is not due to recoils, so the rate of events is independent of their velocity and depends only on the total number density of dark matter particles, thermalised in the vicinity of the detector. In the cases considered here, the effect of gravitational focusing should alter the time variation of the thermalisation depth distribution (see Appendix C). We haven’t taken this effect into account in our simulations, but if we assume that most of the particles thermalised at various depths have equal chances of getting inside the detector, than the time shift of the signal should be $\lesssim 20$ days [29]. However, a typical time delay of the annual modulation that we observe for SIMPs varies from $60$ to $200$ days, so even the maximal correction provided by gravitational focusing cannot improve the fit of the DAMA data sufficiently. Furthermore, the expected problematic diurnal modulation in case of SIMPs is not affected by the gravitational focusing of the Sun, although it can be influenced by that of the Earth. The implication of the latter effect for direct-detection searches was studied in [31]. It appears that the Earth gravitational effect on the diurnal modulation of the total density of dark matter is comparable to the diurnal modulation coming from the rotation of the Earth w.r.t the Galactic frame and, thus, can hardly solve this issue for SIMPs.

5 Conclusion

In this paper, we have shown that the time dependence of the DAMA signal is incompatible with the hypothesis that it comes from SIMPs. Of course, this does not rule out the possibility that SIMPs exist, but they cannot be the source of the signal that DAMA observes. One has then to rely on models with very specific properties tuned to the data, if the signal is due to dark matter.

Soon, this puzzle will have one more chapter as several experiments are under way or planned [32, 33, 34, 35, 36, 37] using the same materials and techniques as DAMA does. First results from COSINE-100 [34] seem to suggest that the DAMA signal is not reproduced, and rule out spin-independent interactions as the cause. KIMS [35] and NAIAD [38] only rule out part of the region of parameter space compatible with DAMA. The data from DM-Ice [32] and ANAIS [36] is consistent with a null-signal, though their statistics is insufficient to constrain the considered parameter region. Note that dark matter must go through the Earth to reach the South Pole, so a potential signal in DM-Ice could not be caused by SIMPs.

As these experiments are rather delicate, one should wait for more statistics and for more independent measurements to close this chapter.

Acknowledgements

We thank Igor P. Ivanov for pointing out the problem studied in this paper. This work was supported by the Fonds de la Recherche Scientifique - FNRS, Belgium, under grant No. 4.4501.15, and M.L. acknowledges the support of FRIA (F.R.S.-FNRS).

Appendix A Solution of the drift-diffusion equation

Let’s consider a one-dimensional drift-diffusion equation in the form

[TABLE]

This equation can be brought to the form

[TABLE]

by the following substitution

[TABLE]

Using this trick one can construct the Green’s function $G^{\prime}(z,\xi,t,\tau)$ for the operator

[TABLE]

given the Green’s function $G(z,\xi,t,\tau)$, which is a solution of the operator

[TABLE]

Similarly to Eq. (13), one has

[TABLE]

Indeed, acting with the operator $L^{\prime}$ on $G^{\prime}$ and performing some simple derivations give

[TABLE]

which is equal to $\delta(z-\xi)\delta(t-\tau)$ for any $z$ and $t$ and, thus, $G^{\prime}$ is the Green’s function for the operator $L^{\prime}$.

The solution of the drift-diffusion equation with the source $f(z,t)$

[TABLE]

is then given by

[TABLE]

where $G(z,\xi,t,\tau)$ is the Green’s function of the simple diffusion equation without the drift term.

The one-dimensional case considered here can be easily generalised to multi-dimensional case.

Appendix B Derivation of the diffusion coefficient

We are going to derive the diffusion coefficient for a gas of particles with mass $M$ spreading through a solid substance, which is comprised of $n$ molecules per unit volume with mass $m$. The latter particles are considered to be significantly lighter than the former and the cross section of their elastic scattering is $\sigma$. The gas is in thermal equilibrium with the surrounding matter at the temperature $T$.

First, let us find the flux of particles $J^{+}$ crossing the plane $z_{0}$ (see Fig. 10) in the direction of $z$. In the context of diffusion, the flux can be understood as the average number of particles $N^{+}$ per unit area, which pass through the given surface in between collisions, divided by the mean free time $\langle\tau\rangle$

[TABLE]

To derive a more useful expression, consider an infinitesimal layer of particles, located at the distance $l$ from the plane $z_{0}$. The number of particles per unit area of the layer $(l,l+dl)$, which move in the direction of z at angles ($\theta,\theta+d\theta)$ is given by

[TABLE]

where the value of $\theta$ ranges from [math] to $\pi/2$, because the particles moving in the opposite direction to $z_{0}$ can never reach it. However, not all of the particles determined by the above expression will manage to pass through $z_{0}$ – some of them will experience a collision before. The distance $S$ that a particle travels before the next encounter is a random variable with an exponential probability density distribution

[TABLE]

where $\lambda=(n\sigma)^{-1}$ is the mean free path. Thus, the probability for particles defined by Eq. (21) to reach the surface at the distance $l/\cos{\theta}$ is

[TABLE]

However, we are interested not in the distance between collisions itself, but in the distance that a particle travels before significantly changing the direction of its velocity. When a heavy particle collides with a much lighter one, as in the case that we consider, it tends to kinematically preserve the direction of its movement (the effect known as the persistence of velocity, see Chap. V in [39]). Suppose, that the probability for a particle to proceed in the same direction after a collision is $q$. Then, the mean path length in one direction is

[TABLE]

where $k$ is a persistence coefficient, which we are going to derive later. Thus, for our purposes we should substitute this increased mean path value in Eq. (22).

Now, gathering Eqs. (21) and (23), we are ready to calculate the numerator in Eq. (20)

[TABLE]

We can expand the density in the vicinity of $z_{0}$ as

[TABLE]

where we have assumed that the density is decreasing in the direction of $z$, so that the net flux of particles through $z_{0}$ in that direction should be positive. Substituting this into Eq. (25) and performing the integration one obtains

[TABLE]

Dividing this by the mean free time

[TABLE]

where $P(v)$ is the Maxwell-Boltzmann distribution for speeds and $\langle v\rangle=\sqrt{8kT/\pi M}$ is the mean speed of gas particles under consideration, we arrive at the expression for $J^{+}$

[TABLE]

The flux $J^{-}$ crossing the plane $z_{0}$ in the opposite direction is given by the same formula, but the sign of the second term is positive due to the fact, that we redefine $l$ to be increasing in the same direction as $z$ and the sign in Eq. (26) changes. Then, the net flux $J$ through the surface $z_{0}$ is

[TABLE]

Using Fick’s first law of diffusion

[TABLE]

we derive the expression for the diffusion coefficient

[TABLE]

Now, let us derive the persistence coefficient $k$. The probability $q$ of a particle to preserve the direction of its velocity after a random collision can be regarded as the value averaged over all collisions of the projection of the velocity after collision $\vec{v_{f}}$ onto the direction of its initial velocity $\vec{v_{i}}$, divided by the absolute value of the latter. For the particles that we consider, $\vec{v_{f}}$ can be derived from energy-momentum conservation

[TABLE]

where $\vec{u}$ is the velocity of a target particle, $v_{\rm rel}=|\vec{u}-\vec{v_{i}}|$ and $\vec{n}$ is a random unit vector. We shall first find the averaged projection over all possible directions of $\vec{u}$ and $\vec{n}$, while keeping the absolute values of $\vec{v}$ and $\vec{u}$ fixed. Clearly, the first term in (33) averages to [math], so that the expectation of the component of velocity of either particle after collision in any direction is equal to the component of $V_{\rm CM}$ in that direction. We should not, though, suppose that all directions of $\vec{u}$ are equally likely, because the probability of collision with any two velocities is proportional to $v_{\rm rel}$. Thus the value of the component of $\vec{v_{f}}$ in the direction of $\vec{v_{i}}$ averaged over all possible kinematic directions is

[TABLE]

where $\omega$ is the angle between $\vec{v_{i}}$ and $\vec{u}$. Then, the probability $q$ can be expressed as a function of speeds

[TABLE]

where $v_{\rm rel}$ varies from $|v_{i}-u|$ to $v_{i}+u$. Since we assume that $M\gg m$, it is reasonable then to assume also that $u>v_{i}$ most of the time. Performing the integration and combining the second and the third terms, we get

[TABLE]

Taking into account again, that $M\gg m$, we can suppose that on average $u\gg v$ and hence

[TABLE]

We see, indeed, that the probability for a massive particle to follow the same direction after collision with a much lighter particle is very close to $1$, so the mean free path length in one direction is extended. Finally, from Eq. (24) and Eq. (37) we can write down the persistence coefficient $k$

[TABLE]

where $\mu$ is the reduced mass. At last, plugging this result back into Eq. (32) one obtains the expression for the diffusion coefficient

[TABLE]

Interestingly enough, there is another way to obtain this expression. One should consider the gas of heavy particles in a force-field, derive its mean (drift) velocity in this field by averaging over all possible directions in a random collision and apply the Einstein relation.

Appendix C Dark-matter velocity and thermalisation depth distributions

We assume that SIMPs follow the same velocity distribution as the rest of dark matter. For the dark-matter velocity distribution in the halo rest frame we take a bounded Maxwellian distribution

[TABLE]

where $\theta$ denotes the Heaviside step-function, $v_{0}$ is the dispersion velocity of the halo and $v_{\rm esc}$ as the Galactic escape velocity at the position of the Sun. We adopt the following values of these parameters $v_{0}=220{\,\rm km}/{\,\rm s}$ [40] and $v_{\rm esc}=544{\,\rm km}/{\,\rm s}$ [41]. Here we omit the normalisation factor of the distribution for simplicity. In the laboratory frame dark matter velocity distribution is different

[TABLE]

where $\vec{v}_{\rm lab}$ is the lab velocity in the Galactic frame. This velocity is a periodical function of time and depends on the location on the globe. For our analysis we calculate $\vec{v}_{\rm lab}$ at the location of LNGS following the procedure described in [18].

Using the relation (1) between the velocity $v$ of the dark-matter particle and the distance $l$ it travels in the ground before thermalisation one can transform the velocity distribution $\omega_{v}(\vec{v},t)$ into the thermalisation depth distribution $\omega_{l}(\vec{l},t)$. Here $\vec{l}$ denotes the vector which connects the point where dark-matter particle enters the ground and the point where it thermalises. The direction of $\vec{l}$ obviously coincides with the direction of $\vec{v}$, so the components of these vectors are simply related as $l_{i}=(l/v)v_{i}$. Following the probability density equality

[TABLE]

where

[TABLE]

is the Jacobian of the transformation from the velocity space to the distance space. Thus, for the thermalisation depth distribution one obtains

[TABLE]

This leads to distributions illustrated in Fig. 11.