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

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.


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 O(100) GeV. Direct detection of the DM particles has been the target of numerous experimental programs [1][2][3]. 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 O(10) − O(100) GeV, while at smaller masses the existing experimental constraints [4] are less severe. This has motivated an increasing interest to the low mass region m DM 1 GeV, see e.g. [5] 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 [6], 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 [6], 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 [7,8], 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 [9] 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 O(1) GeV, wherein the direct detection technique described in [6] 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 v 1 . The double-differential recoil rate per unit detector mass is [10] where m DM and m N are the masses of the DM particle and nucleus, respectively. The local DM density is denoted by ρ 0 = 0.3 GeV/cm 3 and the direction of the recoiling nucleus by the unit vectorq. The squared scattering matrix element |M| 2 is summed and averaged over the inital and final spins. The δ-function imposes the kinematic of the elastic scattering, and the minimum WIMP speed required to excite a nuclear recoil of energy where µ DM,N is the reduced mass of the DM-nucleus system. The angular differential event rate is obtained from (1) by integrating over energy: where E min is the threshold energy for creating a detectable recoil event, and E max is the maximum energy allowed by the event selection of the experiment. If no upper bound is imposed by the detection technique, E 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 max = ∞ 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 (v). In this work we will use the Standard Halo Model, defined as a truncated Maxwellian distribution where v e = 537 km/s is the escape velocity, the circular speed v 0 = 220 km/s and the normalization constant is given by Taking all the above together, the angular differential rate becomes The integration over the recoil energy E and the DM velocity v is affected by the fact that the squared matrix element can in principle depend both on q and v. For a systematical analysis, we consider the non-relativistic effective field theory constructed in [11]. 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 which by construction satisfies v ⊥ · 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 [8,11,12] 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 ∼ 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 asf corresponds to the velocity dependence 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 2 ⊥ , and leads to the transverse Radon transform The angular differential event rate (6) can then be expanded as where A is the mass number of the nucleus and σ 0 = 1/(16πA 2 (m DM + m 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 [8]. 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 [5].

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 [13,14] carried out in Ref. [6]. 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 θ angle was selected as cos −1 (1 − 2u) where u is a uniformly distributed random number between 0 and 1 [15]. 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 [13,[16][17][18] 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. [6], 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 = dΩ(dR/dΩ) directly as a Monte Carlo integral over the solid angle Ω, as explained in more detail below.

A. 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 min = 20 eV) of the velocity distribution f 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 θ = π/4. We notice that for a small DM mass, both functions are strongly peaked towards forward recoil, θ = 0, and hence the behavior of the angular differential rate in the low-mass region will be similar regardless of the v 2 ⊥ -dependence of the squared matrix element |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 As the Earth rotates around its axis, the direction of the DM wind with respect to the labframe 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 = dΩ(dR/dΩ) by a Monte Carlo integral over the solid angle Ω, 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 (θ i , φ i ) on the surface of the unit-sphere, with the corresponding value for the threshold energy E min (θ i , φ i ) obtained from the list. For each (θ i , φ i )-point this procedure yields the corresponding differential event rate dR(θ i , φ i ). The total event rate is then obtained as the sum over the points in the list: where the dependence on time t follows from the time-dependence of the laboratory's motion in the galactic rest frame 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 points = 5000 the result of the sum is within 3% of the result for using the total ∼ 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 DM = 0.3 GeV in the top left figure, is very similar for the v 0 -and v 2 ⊥ -interactions. As the DM mass is increased, the expected modulations become more different, but the amplitude of the signal quickly drops below O(1%), and thus undetectable. However, for the DM mass, 340 MeV m DM 450 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 = a 2 n + b 2 n , where 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 2 ⊥ -interaction in red, as a function of the DM mass. For any value of the DM mass above m DM 340 MeV at least one of the ratios is substantially different to allow separation of the v 0 and v 2 ⊥ interactions, as long as the Fourier-components can be reliably reconstructed from the data. Above m DM 450 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.     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 where R is the average event rate over the time interval ∆t, and R(t) is the event rate as a function of time. Figure 7 shows R RMS as a function of the WIMP mass, for Germanium and Silicon.

B. 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 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 MeV m DM 450 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 min depends on the recoil direction. As established in [6], 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 [6] 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 distribu-     tion 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 MeV m DM 450 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 MeV m DM 350 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 [20] and single electron threshold has already been demonstrated in phonon mediated detectors [21]. 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.
where V z = V ·ŵ. For a Maxwell distribution f M (v) = (2πσ 2 v ) − 3 2 exp(−v 2 /(2σ 2 v ) these are explicitly given as:f The SHM distribution (4) where v e is the escape velocity. Then the Radon transforms are given as: The directional event rate (3) is obtained by integration over energy: where we now denote V z = V ·q and m = 2µ 2 DM,N /m N . For the transverse Radon transform the integral over energy reads: