Velocity Dependent Dark Matter Interactions in Single-Electron Resolution Semiconductor Detectors with Directional Sensitivity

TL;DR

This paper explores how velocity and momentum dependence of dark matter interactions influence signals in single-electron resolution semiconductor detectors, revealing daily modulation effects sensitive to dark matter properties.

Contribution

It provides a detailed analysis of velocity and recoil momentum dependence in dark matter interactions and identifies conditions under which daily modulation signals can be observed in semiconductor detectors.

Findings

Dark matter interactions cause daily modulation signals in specific mass ranges.

The structure of modulation is sensitive to velocity and momentum dependence.

Optimal mass ranges for observing these effects are identified.

Abstract

We investigate the velocity and recoil momentum dependence of dark matter interactions with ordinary matter. In particular we focus on the single-electron resolution semiconductor detectors, which allow experimental assessment of sub-GeV dark matter masses. We find that, within a specific mass range depending on the detector material, the dark matter interactions result in a signal characterized by daily modulation. Furthermore, we find that the detailed structure of this modulation is sensitive to the velocity and momentum dependence of dark matter interactions. We identify the optimal mass range for the prevalence of these effects.

Figures

Equations

$$\frac{d^{2}R}{dEd{\Omega }_q}=\frac{{\rho }_0}{2\pi m_{\mathrm{DM}}}\frac{|\mathcal{M}{|}^2}{32\pi m_{\mathrm{N}}^2{m}_{\mathrm{DM}}^2}\delta (\mathbf{v}\cdot \hat{\mathbf{q}}-v_{\min }),$$

$$v_{\min }=\sqrt{\frac{m_{\mathrm{N}}E}{2{\mu }_{\mathrm{DM},\mathrm{N}}^2}},$$

$$\frac{dR}{d{\Omega }_q}={\int }_{E_{\min }}^{E_{\max }}\frac{d^{2}R}{dEd{\Omega }_q}dE,$$

$$f_{\mathrm{SHM}}(v)=N_{\mathrm{e}}^{-1}{f}_M(v)\Theta (v_e-v),$$

$$N_{\mathrm{e}}=\mathrm{erf}\left(\frac{v_{\mathrm{e}}}{\sqrt{2{\sigma }_v^2}}\right)-\sqrt{\frac{2}{\pi {\sigma }_v^2}}{v}_{\mathrm{e}}{e}^{-\frac{v_{\mathrm{e}}^2}{2{\sigma }_v^2}}.$$

$$\frac{dR}{d{\Omega }_q}=\frac{{\rho }_0}{2\pi m_{\mathrm{DM}}}\frac{1}{32\pi m_{\mathrm{N}}^2{m}_{\mathrm{DM}}^2}\underset{E_{\min }}{\overset{E_{\max }}{\int }}dE\int d^{3}v|\mathcal{M}{|}^2{f}_{\mathrm{SHM}}(v)\delta (\mathbf{v}\cdot \hat{\mathbf{q}}-v_{\min }).$$

$$v_{\perp }^2=v^2-\frac{q^2}{4{\mu }_{\mathrm{DM},\mathrm{N}}^2},$$

$$|\mathcal{M}{|}^2=a_{1}1+a_2{q}^2+a_3{q}^4+b_1{v}_{\perp }^2+b_2{q}^2{v}_{\perp }^2+b_3{q}^4{v}_{\perp }^2+\cdots$$

$$\hat{f}(v_{\min },\hat{q})=\int d^{3}vf(v)\delta (\mathbf{v}\cdot \hat{\mathbf{q}}-v_{\min }),$$

$${\hat{f}}^{\mathrm{T}}(v_{\min },\hat{q})=\int d^{3}vf(v)v_{\perp }^2\delta (\mathbf{v}\cdot \hat{\mathbf{q}}-v_{\min }).$$

$$\frac{dR}{d{\Omega }_q}=\frac{{\rho }_0}{4\pi m_{\mathrm{DM}}}\frac{{\sigma }_0{A}^2}{{\mu }_{\mathrm{DM},\mathrm{N}}^2}\underset{E_{\min }}{\overset{E_{\max }}{\int }}dE\left((a_1+a_2{q}^2+\dots ){\hat{f}}_{\mathrm{SHM}}(v_{\min },\hat{q})+(b_1+b_2{q}^2+\dots ){\hat{f}}_{\mathrm{SHM}}^{\mathrm{T}}(v_{\min },\hat{q})\right),$$

$$R(t)=\frac{4\pi }{N_{\mathrm{points}}}\sum _{i=1}^{N_{\mathrm{points}}}dR({\theta }_i,{\varphi }_i,t),$$

$$a_n=\underset{0}{\overset{1}{\int }}R(t)\cos \left(2\pi nt\right)dt,b_n=\underset{0}{\overset{1}{\int }}R(t)\sin \left(2\pi nt\right)dt,$$

$$R_{\mathrm{RMS}}=\sqrt{\frac{1}{⟨R{⟩}^2\Delta t}{\int }_{\Delta t}(R(t)-⟨R⟩{)}^{2}dt},$$

$$\hat{f}(w,\hat{\mathbf{w}})=\int \delta (w-\hat{\mathbf{w}}\cdot \mathbf{v})f(\mathbf{v})d^{3}v,$$

$${\hat{f}}^T(w,\hat{\mathbf{w}})=\int \delta (w-\hat{\mathbf{w}}\cdot \mathbf{v})({\mathbf{v}}_{\perp }{)}^{2}f(\mathbf{v})d^{3}v,$$

$$\hat{f}(w,\hat{\mathbf{w}})$$

$${\hat{f}}^T(w,\hat{\mathbf{w}})$$

$$\hat{f}(w,\hat{\mathbf{w}})$$

$${\hat{f}}^T(w,\hat{\mathbf{w}})$$

$${\hat{f}}_M(w,\hat{\mathbf{w}})$$

$${\hat{f}}_M^T(w,\hat{\mathbf{w}})$$

$${\hat{f}}_{\mathrm{SHM}}(w,\hat{\mathbf{w}})$$

$${\hat{f}}_{\mathrm{SHM}}^T(w,\hat{\mathbf{w}})$$

$$\underset{E_{\min }}{\overset{\infty }{\int }}$$

Full text

††preprint: HIP-2019-8/TH

I Introduction

Cosmological and astrophysical observations provide overwhelming evidence for the existence of dark matter (DM) consisted of particles beyond the Standard Model (SM). The well established paradigm is a Weakly Interacting Massive Particle (WIMP), with an electroweak scale mass $\mathcal{O}(100)$ GeV. Direct detection of the DM particles has been the target of numerous experimental programs Akerib:2017kat ; Cui:2017nnn ; Aprile:2018dbl . So far these experiments have not led to a consistent discovery, but have provided solid constraints limiting the strength of the interactions between dark and ordinary matter. The current direct detection experiments are most effective around the typical WIMP mass range of $\mathcal{O}(10)-\mathcal{O}(100)$ GeV, while at smaller masses the existing experimental constraints Agnese:2017njq are less severe. This has motivated an increasing interest to the low mass region $m_{\rm DM}\lesssim 1$ GeV, see e.g. Battaglieri:2017aum and the references therein. Since the coherent elastic scattering of solar neutrinos produce signals that mimic those expected from low mass DM interactions, these experiment will eventually become background limited and need to develop methods to mitigate this irreducible background.

A promising method for direct detection of low mass DM particles was presented in Kadribasic:2017obi , based on single-electron resolution semiconductor detectors. Assuming a direct correlation between ionization and defect creation thresholds in semiconductors, it was noted that due to the anisotropic structure of the semiconductor crystals, the ionization threshold can be sensitive to the recoil direction. Hence this technique allows for a detection of a daily modulation signal, due to the rotation of the earth with respect to the direction of the DM wind.

The purpose of this paper is to extend the analysis presented in Kadribasic:2017obi , to cover the variety of non-relativistic operators describing the DM-SM scattering. The directional detection in the context of non-relativistic effective theory of DM-SM scattering has been discussed in Catena:2015vpa ; Kavanagh:2015jma , where the angular recoil distributions expected for the various effective operators have been described. Recording the full angular recoil spectrum could thus be used to identify the operator characterizing the DM-SM scattering, and would reveal valuable information about the underlying theory of DM. See Mayet:2016zxu for an overview of the prospects in directional DM detection.

Our method, however, does not rely on recording the angular differential recoil distribution, but rather in observing the daily modulation in the total integrated event rate, the origin of which is in the directional sensitivity of the threshold energy. Therefore, the question becomes, to which extent is the modulation signal sensitive to the type of the effective DM-SM scattering operator? Qualitatively, the following behavior is to be expected: The amplitude of the daily modulation signal grows towards lighter DM mass, as the threshold energy for creating the electron-hole pair becomes more significant in comparison to the average kinetic energy available for the recoil. On the other hand, as the DM mass is decreased, the angular recoil spectrum becomes peaked at forward direction for all effective operators. This effect is due to the recoil kinematics i.e. if the recoil angle is large, only a small fraction of the kinetic energy of the incoming DM particle is available for the recoil energy of the nucleus.

In the following we will quantify the above assertions, and provide analytic formulas for the relevant event rates. Based on our results, we establish the following general picture i.e. in the small mass region where the daily modulation amplitude is readily observable, the recoil dynamics are largely insensitive to the effective DM-SM scattering operator. As the DM mass increases, the different recoil dynamics cause an increasing difference in the angular distribution. However, this effect becomes masked by the fact that the amplitude of the daily modulation signal quickly decreases as the DM mass increases, thus the effect becomes less observable. We will thus arrive to the conclusion, that depending on the detector material, there is a range of DM masses below $\mathcal{O}(1)$ GeV, wherein the direct detection technique described in Kadribasic:2017obi is valid. At the lower edge of this mass range, the shape of the modulation signal is practically blind to the details of the underlying theory of DM-SM scattering, but for increasing masses the shape of the modulation signal strongly depends on the velocity dependence of the underlying DM-SM interactions. The most important parameter for the efficiency of this method is the DM mass.

The paper is organized as follows: in section II we introduce the kinematics of the scattering event and the non-relativistic operators. In section III we describe the directional sensitivity of our detector concept and the resulting daily modulation signal, and we conclude in section IV. The analytic formulas for the angular event rates are given in appendix A.

II Basic formulas

Consider scattering of a dark matter particle with a nucleus. Denote the DM velocity by ${\bf{v}}$111Throughout this paper we will denote three-vectors, such as ${\bf{v}}$, with a boldface font and their amplitudes with italic, as $v$.. The double-differential recoil rate per unit detector mass is Gondolo:2002np

[TABLE]

where $m_{\rm{DM}}$ and $m_{\rm N}$ are the masses of the DM particle and nucleus, respectively. The local DM density is denoted by $\rho_{0}=0.3\ {\rm GeV}/{\rm cm}^{3}$ and the direction of the recoiling nucleus by the unit vector $\hat{{\bf{q}}}$. The squared scattering matrix element $|{\cal{M}}|^{2}$ is summed and averaged over the inital and final spins. The $\delta$-function imposes the kinematic of the elastic scattering, and the minimum WIMP speed required to excite a nuclear recoil of energy $E=q^{2}/2m_{\rm N}$ is

[TABLE]

where $\mu_{{\rm{DM,N}}}$ is the reduced mass of the DM-nucleus system. The angular differential event rate is obtained from (1) by integrating over energy:

[TABLE]

where $E_{\rm min}$ is the threshold energy for creating a detectable recoil event, and $E_{\rm max}$ is the maximum energy allowed by the event selection of the experiment. If no upper bound is imposed by the detection technique, $E_{\rm max}$ can be taken to infinity, as the convergence of the integral is ensured by the integrability of the DM velocity distribution, to be discussed below. In this paper we will take $E_{\rm max}=\infty$ unless otherwise noted.

To calculate the observable directional event rate in a detector on earth, the distribution of DM velocities in the galactic halo must be taken into account: the rate in Eq.(1) must be integrated over all DM velocities weighted by the distribution $f({\bf{v}})$. In this work we will use the Standard Halo Model, defined as a truncated Maxwellian distribution

[TABLE]

where $v_{\rm e}=537\ {\rm km}/{\rm s}$ is the escape velocity, $f_{M}(v)=(2\pi\sigma_{v}^{2})^{-3/2}\exp(-v^{2}/2\sigma_{v}^{2})$ is the Maxwellian distribution with a standard deviation $\sigma_{v}=v_{0}/\sqrt{2}$, the circular speed $v_{0}=220\ {\rm km}/{\rm s}$ and the normalization constant is given by

[TABLE]

Taking all the above together, the angular differential rate becomes

[TABLE]

The integration over the recoil energy $E$ and the DM velocity ${\bf{v}}$ is affected by the fact that the squared matrix element can in principle depend both on $q$ and ${\bf{v}}$. For a systematical analysis, we consider the non-relativistic effective field theory constructed in Fitzpatrick:2012ix . The effective field theory operator basis is constructed by imposing the requirement of Hermiticity together with invariance under Galilean transformations and time reversal. In particular, because of Hermiticity, the velocity dependence of these operators is only through the combination

[TABLE]

which by construction satisfies ${\bf{v}}_{\perp}\cdot{\bf{q}}=0$.

For our purposes it is sufficient to categorize different interactions in terms of the velocity and energy dependence they imply for the square of the averaged matrix element appearing in Eq. (6). The possible dependences are Kavanagh:2015jma ; Fitzpatrick:2012ix ; DelNobile:2018dfg

[TABLE]

where the ellipsis stands for operators of higher order in $q^{2}$ and $a_{i},b_{i}$ are coefficients with mass dimension $-2(i-1)$. In addition to these, it is interesting to consider effects from long rage interactions mediated by some light field. These will lead to behavior $\sim q^{-4}$.

Hence, in order to probe the full range of different behaviors due to different interactions, we only need to compute two different integrals over the velocity distribution. The Radon transform, defined as

[TABLE]

corresponds to the velocity dependence ${\cal{O}}(v^{0})$ of the matrix element. The only other possibility, then, is that the squared matrix element is proportional to the square of the perpendicular velocity $v_{\perp}^{2}$, and leads to the transverse Radon transform

[TABLE]

The angular differential event rate (6) can then be expanded as

[TABLE]

where $A$ is the mass number of the nucleus and $\sigma_{0}=1/(16\pi A^{2}(m_{\rm DM}+m_{\rm N})^{2})$ is a reference DM-nucleon cross section. For the SHM the integral over energy can be performed analytically, and the necessary explicit formulas are provided in the appendix A.

Notice that the overall normalization of the terms in the expansion (11) must include the corresponding nuclear matrix elements Kavanagh:2015jma . In this work our goal is not to determine these absolute normalizations, but to determine the shape of the resulting observable signal which our detector concept would measure given enough exposure, and whether the shape of the signal, as a function of time, is sensitive to the structure of the underlying operators. Therefore we absorb the normalization of the operators in the coefficients $a_{i},b_{i}$. We also neglect the nuclear form-factors which suppress high energy recoils. In the low-mass region we are considering, the recoil energies are small and the form factors are very close to one. We provide a compilation of necessary general formulas for the event rates, which complement existing literature and are expected to be useful for similar studies within different detector concepts currently under active investigation Battaglieri:2017aum .

III Directional Energy Threshold

In semiconductor materials the threshold energy for defect creation is a function of the recoil direction. The representation of this effect for Germanium and Silicon is shown in figure 1. To obtain this figure we have generated a sample of 84936 randomly sampled directions in Germanium and 24155 directions in Silicon, with the corresponding energy thresholds, utilizing the data from the molecular dynamics simulations Holmstrom:2008 ; Nordlund:2006 carried out in Ref. Kadribasic:2017obi . Briefly, Ge and Si atom recoils were simulated in randomly generated directions in three dimensions. The directions were selected to give a uniform distribution over solid angle, i.e. the $\theta$ angle was selected as $\cos^{-1}(1-2u)$ where $u$ is a uniformly distributed random number between 0 and 1 Press:1992zz . For each direction, the recoil energy was increased from 4 eV in 1 eV increments until a stable defect was produced. Time-dependent density functional theory calculations Lim:2016 ; Horsfield:2016 ; Holmstrom:2008 ; Holmstrom:2010 showed that also the ionization has a strong dependence on crystal directions. Unfortunately these calculations are too demanding computationally to obtain a full threshold map, and hence we continue to work with the inference that the ionization energy threshold correlates with the defect production threshold.

Consequently, the event rates obtained by integrating the Radon transforms (9) and (10) over energy become functions of the recoil direction. Contrary to Ref. Kadribasic:2017obi , in the current work we did not average the threshold energy surface over an angular interval. Instead, we use the list of randomly sampled directions with the corresponding energy thresholds to compute the event rate $R=\int d\Omega(dR/d\Omega)$ directly as a Monte Carlo integral over the solid angle $\Omega$, as explained in more detail below.

III.1 Purely velocity dependent interactions

To understand the general behavior of the angular differential rate as a function of the DM mass, we begin by showing the integrated Radon transforms (for $E_{\rm min}=20$ eV) of the velocity distribution $f_{\rm SHM}$, in figure 2 for various values of the WIMP mass. For the purpose of illustration, the functions have been arbitrarily normalized so that they match at the point $\theta=\pi/4$. We notice that for a small DM mass, both functions are strongly peaked towards forward recoil, $\theta=0$, and hence the behavior of the angular differential rate in the low-mass region will be similar regardless of the $v_{\perp}^{2}$-dependence of the squared matrix element $|{\cal M}|^{2}$. For larger values of the DM mass both distributions become broader, and the transverse Radon transform develops a maximum at some non-zero recoil angle.

Then, to demonstrate the effect of the direction-dependent energy threshold, we show the angular event rate for the $v^{0}$- (the Radon transform) and $v_{\perp}^{2}$- (the transverse Radon transform) interactions in figure 3 for various values of the WIMP mass. These figures are obtained by integrating the Radon transforms over energy, with $E_{\rm min}$ in each direction given by the data shown in figure 1. For these calculations, we have assumed the SHM velocity distribution and a Germanium detector on the SNOLAB site $(46.4719^{\circ}{\rm N},81.1868^{\circ}{\rm W})$ on September 6, 2015 at 18:00. The event rates correspond to the Spin-independent DM-nucleon cross section $a_{1}\sigma_{0}=10^{-39}\ {\rm cm}^{2}$ in the $v^{0}$-case, and $b_{1}\sigma_{0}=10^{-33}\ {\rm cm}^{2}$ for the $v_{\perp}^{2}$-interaction.

We notice that the distributions for the small DM masses, shown in the top row of the figure, are basically indistinguishable by eye, and centered towards the average direction of the DM wind. As the DM mass is increased, the angular recoil distribution becomes wider. Eventually, for large enough $m_{\rm DM}$, the information of the directionality of the energy threshold becomes practically undetectable, as is evident from the figures on the bottom row. Towards large $m_{\rm DM}$, the off-zero maximum of the transverse Radon transform manifests as the ring-like feature around the direction of the DM wind, visible in the bottom right figure.

As the Earth rotates around its axis, the direction of the DM wind with respect to the lab-frame modulates. Due to the directional dependence of the event rate shown above, this results in a diurnal modulation of the integrated event rate. Figure 4 shows the diurnal modulation in the event rate for various values of the WIMP mass during September 6, 2015. We compute the event rate $R=\int d\Omega(dR/d\Omega)$ by a Monte Carlo integral over the solid angle $\Omega$, utilizing the list of 84936 randomly sampled directions with the corresponding energy thresholds. Energy integrals of the Radon transforms are evaluated for each sampled point $(\theta_{i},\phi_{i})$ on the surface of the unit-sphere, with the corresponding value for the threshold energy $E_{\rm min}(\theta_{i},\phi_{i})$ obtained from the list. For each $(\theta_{i},\phi_{i})$-point this procedure yields the corresponding differential event rate $dR(\theta_{i},\phi_{i})$. The total event rate is then obtained as the sum over the points in the list:

[TABLE]

where the dependence on time $t$ follows from the time-dependence of the laboratory’s motion in the galactic rest frame ${\bf{V}}(t)$, as explained in the appendix A. We have checked that the number of points is sufficient for an accurate integral: Already for $N_{\rm points}=5000$ the result of the sum is within 3% of the result for using the total $\sim 85000$ points in the list.

As expected from the discussion above, the shape of the diurnal modulation signal for the smallest DM mass displayed here, $m_{\rm DM}=0.3$ GeV in the top left figure, is very similar for the $v^{0}$- and $v_{\perp}^{2}$-interactions. As the DM mass is increased, the expected modulations become more different, but the amplitude of the signal quickly drops below $\mathcal{O}(1\%)$, and thus undetectable. However, for the DM mass, $340\ {\rm MeV}\lesssim m_{\rm DM}\lesssim 450\ {\rm MeV}$, the shape of the daily modulation signal can be used to determine the velocity-dependence of the scattering amplitude.

To assess the feasibility of velocity-dependence detection from the shape of the daily modulation signal, we analyze the Fourier-components of the daily event rates, $C_{n}=\sqrt{a_{n}^{2}+b_{n}^{2}}$, where

[TABLE]

where $t$ is time in units of day. We show the ratios $C_{1}/C_{2}$, $C_{3}/C_{2}$ and $C_{4}/C_{2}$ in figure 5, for the $v^{0}$-interaction in black and $v_{\perp}^{2}$-interaction in red, as a function of the DM mass. For any value of the DM mass above $m_{\rm DM}\gtrsim 340\ {\rm MeV}$ at least one of the ratios is substantially different to allow separation of the $v^{0}$ and $v_{\perp}^{2}$ interactions, as long as the Fourier-components can be reliably reconstructed from the data. Above $m_{\rm DM}\gtrsim 450\ {\rm MeV}$ the amplitude of the daily modulation rate drops below 1%, and the reconstruction of the Fourier-components becomes prohibiting in the required scale of the experiment. Figure 5 also shows the ratios of the Fourier components for a Silicon detector, where the $v^{0}$-interaction is shown by the gray dashed line, and the $v_{\perp}^{2}$-interaction by the purple dashed line. Due to the smaller atomic mass of Silicon compared to Germanium, the Si sensitive region falls at the lower values of the DM mass. We identify the range of $250\ {\rm MeV}\lesssim m_{\rm DM}\lesssim 350\ {\rm MeV}$ as the region where the velocity-dependence of the operator can be identified in Silicon.

To understand how the shape of the daily modulation signal arises, we show in figure 6 the event rate for the $v^{0}$-interaction as a function of the recoil direction, for $m_{DM}=0.3$ GeV, corresponding to the top left panel of figure 4, at the moments of minimum and maximum event rates (at 04:00, 10:00, 16:00, 22:00 hours), assuming $a_{1}\sigma_{0}=10^{-39}\ {\rm cm}^{2}$. Comparing to figure 1, we see that the maximum event rates (corresponding to 10:00 and 22:00), shown on the right column of figure 6, occur when the direction of the DM wind, shown by the blue dot in figure 6, coincides with the low threshold energy directions that appear as the dark spots in figure 1. Respectively, the minima of the event rate (at 4:00 and 16:00, shown on the left) occur when the direction of the DM wind is maximally far away from the low threshold regions. The blue curve in the figure shows the path of the direction of the DM wind on the unit sphere during the 24 hour period.

We define the normalized RMS daily modulation by

[TABLE]

where $\langle R\rangle$ is the average event rate over the time interval $\Delta t$, and $R(t)$ is the event rate as a function of time. Figure 7 shows $R_{\rm RMS}$ as a function of the WIMP mass, for Germanium and Silicon.

III.2 Energy dependent interactions

To demonstrate the behavior of the $q^{2}$-dependent scattering operators, we will focus here on the leading term in the $q^{2}$-expansion, the $q^{2}$, and the long-range force effective operator, $q^{-4}$.

Figure 8 shows the integrated Radon transform (with $E_{\rm min}=20$ eV) as a function of the recoil direction for the operators $1$ (black), $q^{2}$ (red) and $q^{-4}$ (blue), for various values of the DM mass. Again we notice that the functions become similar to each other for small values of the DM mass.

This is also apparent in the daily rates, shown in figure 9, for $m_{DM}=0.3$ GeV and $m_{DM}=0.4$ GeV wherein, the normalized event rates are nearly equal for the 0.3 GeV particle, but begin to deviate for larger values of the DM mass. The RMS modulation as a function of the DM mass is consistent with this observation, as shown in figure 10.

Finally, in figure 11 we show the same ratios of the Fourier-components as was shown in figure 5 but here we also include the long-range interaction $q^{-4}$ and the $q^{2}$-interaction, presented by the purple dashed and the blue dotted lines accordingly. We conclude that within the range $340\ {\rm MeV}\lesssim m_{\rm DM}\lesssim 450\ {\rm MeV}$ identified above, the long-range interaction can also be identified based on the ratios of the Fourier components. The case of the $q^{2}$-interaction shown in blue is more subtle, as the ratios of the Fourier-components for this interaction resemble those for the $v^{0}$-operator. However, within the most promising mass-window, these two can be additionally separated by the use of $C_{1}/C_{2}$-ratio.

IV Summary and outlook

We have considered dark matter scattering in the single-electron resolution ionization detectors wherein the quantum of electronic excitation $E_{\rm min}$ depends on the recoil direction. As established in Kadribasic:2017obi , the signal of dark matter scattering in this case is detectable via the observation of diurnal modulation in expected event rates. We have extended the analysis of Kadribasic:2017obi to cover the possible velocity and energy dependencies of dark matter scattering on ordinary matter as implied by the general low energy effective theory of dark matter.

We carried out the analysis using the Standard Halo Model for the DM velocity distribution and our main finding is that for a given detector material there is a range of sub-GeV masses, wherein the modulation signal is strong. Furthermore, near the lower boundary of this mass range the shape of the modulation signal is practically independent on the nature of the underlying DM interaction. However, above this lower boundary the signal develops a strong dependence on the DM velocity and scattering energy, allowing for discrimination of different classes of interaction operators. At higher DM masses, towards one GeV, the overall amplitude of the modulation signal decreases and becomes less discernable.

In this study we have focused on germanium as the detector material, for which we identified the mass interval $340\ {\rm MeV}\lesssim m_{\rm DM}\lesssim 450\ {\rm MeV}$ as the most promising window where the type of the DM-SM scattering operator can be identified from the shape of the daily modulation signal. We have also performed a preliminary study on silicon, where we find qualitatively similar behavior, and identify the separation window as $250\ {\rm MeV}\lesssim m_{\rm DM}\lesssim 350\ {\rm MeV}$. With a selection of detector materials, it could thus be feasible to cover a larger range of DM masses, with multiple experiments having partly overlapping regions of sensitivity. The exploration of the directional dependence of the ionization energy threshold in a variety of materials is therefore strongly motivated.

There is a rising interest in the dark matter search community to develop very low threshold detectors  contact_free and single electron threshold has already been demonstrated in phonon mediated detectors  CDMS_singleE . A careful calibration of common semiconductors for dark matter detection down to the single electron-hole excitation level is necessary in order to interpret their results. Using mono-energetic neutron beams, such efforts are currently ongoing in various facilities and the results presented in this work can be verified in those experiments.

Appendix A Radon transforms and their energy integrals

Consider the Radon transform of a function $f$, defined as

[TABLE]

and the transverse Radon transform

[TABLE]

where ${\bf{v}}_{\perp}={\bf{v}}-({\bf{v}}\cdot\hat{{\bf{w}}})\hat{{\bf{w}}}$.

Choosing the $z$-axis parallel to the unit vector $\hat{{\bf{w}}}$ and assuming isotropic velocity distribution $f({\bf{v}})=f(v)$, these can be expressed as integral over the amplitude $v$ only:

[TABLE]

The motion with respect to the galactic rest frame is taken into account via the coordinate transformation ${\bf{v}}\rightarrow{\bf{v}}-{\bf{V}}$, where ${\bf{V}}$ is the velocity of the lab-frame with respect to the galactic rest frame. We follow the parametrization given in Bozorgnia:2011tk ; Mayet:2016zxu for ${\bf{V}}(t)$. In this case we have

[TABLE]

where $V_{z}={\bf{V}}\cdot\hat{{\bf{w}}}$. For a Maxwell distribution $f_{M}(v)=(2\pi\sigma_{v}^{2})^{-\frac{3}{2}}\exp(-v^{2}/(2\sigma_{v}^{2})$ these are explicitly given as:

[TABLE]

The SHM distribution (4) is $f_{\rm{SHM}}(v)=N_{\rm{e}}^{-1}f_{M}(v)\Theta(v_{e}-v)$, where $v_{e}$ is the escape velocity. Then the Radon transforms are given as:

[TABLE]

The directional event rate (3) is obtained by integration over energy:

[TABLE]

where we now denote $V_{z}={\bf{V}}\cdot\hat{{\bf{q}}}$ and $m=2\mu_{\rm DM,N}^{2}/m_{\rm N}$. For the transverse Radon transform the integral over energy reads:

[TABLE]

For the $q^{2}$ operator, the Radon transform (23) must be multiplied by $q^{2}$ before taking the integral over energy. The result is

[TABLE]

For the integral $\int\limits_{E_{\rm min}}^{\infty}q^{-4}\hat{f}_{\rm SHM}(v_{\rm{min}},\hat{{\bf{q}}})dE$ we find no analytic expression, and therefore perform the integral over energy numerically in this case.