Impact of uncertainties in the halo velocity profile on direct detection of sub-GeV dark matter

We use the state-of-the-art high-resolution cosmological simulations by IllustrisTNG to derive the velocity distribution and local density of dark matter in galaxies like our Milky Way and find a substantial spread in both quantities. Next we use our findings to examine the sensitivity to the dark matter velocity profile of underground searches using electron scattering in germanium and silicon targets. We find that sub-GeV dark matter search is strongly affected by these uncertainties, unlike nuclear recoil searches for heavier dark matter, especially in multiple electron-hole modes, for which the sensitivity to the scattering cross-section is also weaker. Therefore, by improving the sensitivity to lower ionization thresholds not only projected sensitivities will be boosted but also the dependence on the astrophysical uncertainties will become significantly reduced.


Introduction
Our limited understanding of the dark matter (DM) halo structure, both local and Galactic, constitutes one of the largest sources of uncertainty in direct and indirect DM searches. As astrophysical observations present many challenges in determining the mass, density, shape and velocity profile of dark matter (see, e.g., [1][2][3], where the latter is based on the recent results from the Gaia satellite), the method of choice for studying DM halo properties has often been based on numerical simulations of large scale structures. These, however, until recently suffered from lacking to account for the baryonic component, or an oversimplified implementation of its role. Only very recently progress in hydrodynamical simulations with state-of-the-art baryonic physics effects included allowed for sufficient precision to accurately simulate Milky Way (MW)-like galaxies.
In this paper we use the cosmological simulations of IllustrisTNG [4][5][6] to infer local properties of the DM halo structure of MW-like galaxies and subsequently apply our findings to examine the impact of the estimated astrophysical uncertainties on direct DM searches. Until recently significant effort was devoted to studying the implications of different assumptions about the local DM density and the velocity profile for the direct detection searches of heavy -of a few GeV, or more -DM based on nuclear recoils (see, e.g., a detailed discussion in [7]). Even though our results can be applied to this case as well, in this work we focus on a previously basically unexplored impact on direct detection of sub-GeV dark matter via electron scattering (see [8] for earlier work). As we shall see, in this detection technique DM halo properties, especially its velocity distribution, start to play a significant role.
This paper consists of two main parts: in section 2 we construct and discuss several DM profiles inferred from the Illustris simulation, while in section 3 we apply our findings to study their effect on electron recoil searches of DM with sub-GeV mass in semiconductor targets. Section 4 provides discussion of the results and conclusions.

Dark matter velocity distributions from IllustrisTNG simulation
As outlined in the introduction, dark matter velocity distribution is a crucial ingredient for estimating direct dark matter detection in the Milky Way. Most studies assume that the dark matter halo model follows the Standard Halo Model (SHM) [43], either for simplicity or to provide a platform for comparison among different experiments. The SHM is an isothermal, spherical dark halo with an isotropic Maxwell-Boltzmann (MB) dark matter velocity distribution with the dispersion velocity v 0 (taken to be the local circular speed, hereinafter referred to as peak speed ) in the Galactic rest frame and truncated at the escape speed from the Galaxy, where N is a normalization factor. The local circular speed is usually assumed to be v 0 = 220 km/s and the commonly adopted escape speed is v esc = 544 km/s. A departure from this commonly assumed model changes the predicted event rate in direct detection experiments. In fact, the true dark matter distribution could be different from what is normally assumed. One way to access the information on the properties of the DM halo is to examine the dark matter distribution of the MW analogues in the cosmological simulations. For example, the DM-only simulations of structure formation have shown that the velocity distributions have substantial deviations from a Maxwellian shape [44][45][46][47]. Nevertheless, in recent years by using more realistic high resolution cosmological simulations, which include baryonic physics effects, it was shown that the dark matter velocity distributions actually agree rather well with the SHM. Yet, the local circular speed of the Maxwellian distribution and the escape velocity may well be different [48][49][50][51][52]. In this section, we focus on investigation of the MW-like galaxies from the IllustrisTNG project. We first introduce the general properties of the used simulation and describe the selection criteria used to identify the MW analogues. We then present the resulting velocity distribution functions for some selected halos.

Selection of MW-like galaxies
The IllustrisTNG project is a publicly available suite of state-of-the-art magneto-hydrodynamic cosmological simulations [4][5][6], an extension to the Illustris simulations [53][54][55]. The -2 -JHEP07(2020)081 cosmological parameters were chosen to be those derived from the analysis of the Planck collaboration 2015, XIII [56], namely Ω m = 0.3089, Ω b = 0.0486, h = 0.677, σ 8 = 0.8159 and n S = 0.9667. IllustrisTNG is evolved using the moving mesh hydrodynamics code AREPO [6] and it contains a number of improvements of sub-grid physics models (e.g. models for stellar and AGN feedback, as well as black hole growth [57]) with respect to the original Illustris project. These simulations are able to reproduce many key relations in the observed galaxies, such as the stellar mass function and the size distribution at low and high redshift, the fraction of dark matter within galaxies at z = 0. For more details, see [5,[58][59][60][61][62]. Furthermore, as it was shown in [60], the MW-mass galaxies of the Illus-trisTNG reproduce adequately the observed MW rotation curve, this suggests that the matter distribution in the simulated MW-mass galaxies is realistic and the baryonic effects in the dark matter are properly modelled. That makes them of particular relevance for our study.
In this paper we make use of the highest (to date) resolution version of the IllustrisTNG project, which is TNG100. From the suite of TNG100 simulations we select candidate MW analogues in a (75/h ≈ 110.7 Mpc) 3 box. The baryonic and dark matter mass resolutions are m baryon = 1.4 × 10 6 M and m DM = 7.5 × 10 6 M , respectively.
In order to identify simulated galaxies which satisfy MW observational constraints we adopt the following criteria: • we identify all galaxies in the mass range 1 5 × 10 11 M 200 /M 2 × 10 12 , where M 200 is defined as the mass enclosed within the sphere that contains a mean density 200 times the critical density. The MW mass is usually estimated to be within the assumed mass range, for the compilation of recent MW mass estimates see e.g. [64]; • within 2456 galaxies identified in the previous step we select the ones that have the stellar mass within the following mass range: 2.5 × 10 10 M * /M 5 × 10 10 . The mass range is slightly lower than that usually inferred from MW, see e.g. [65]. This is related to the fact that the effects of black hole feedback are too dominant in TNG simulation, see for more details [63]. We also look at the gaseous mass to stellar mass fraction that we choose to be in the following range: 0.03 M gas /M * 1, as not to have very poor/reach gas galaxies. At this step we end up with 340 galaxies; • among those 340 galaxies we identify galaxies that have a substantial stellar disc component. In order to do so we follow the approach of ref. [66]. We then select only galaxies that contain more than 50% of stars with the circularity parameter > 0.7, i.e., those stars that have significant (positive) rotational support. Here we also require that selected galaxies do not have any significant merger after redshift z = 0.68 and we end up with 164 galaxies.
For each galaxy in the sample, we align the coordinate system with the principal axes. We then select DM particles in torus, along the plane in the stellar disc, around the solar JHEP07(2020)081 radius, R = 8 kpc, with both radial and vertical extents ±2 kpc. This region contains a total of ∼ 700-900 dark matter "particles".

Dark matter velocity distributions
To derive the dark matter distribution we compute the average speed distribution of DM particles in the torus defined above. The speed distribution, f (v), is normalised to unity as dvf (v) = 1, and is related to the velocity distribution,f (v), by Here dΩ v is an infinitesimal solid angle around the direction v. We then use MB distribution introduced in eq. (2.1) to perform a fit to the resultant DM velocity distribution. We find that only ∼ 20% of galaxies in our sample can be well described by the MB velocity distribution, i.e. have the resulting χ 2 ≤ 1, and in our further analysis we restrict ourselves to this choice. While there are arguments that the actual DM velocity distribution function of the MW has a more complicated form (see, e.g. [3,[67][68][69]), we believe that our choice is still sufficient to examine our main task which is to examine the effect of astrophysical uncertainties on projected sensitivities of detecting sub-GeV dark matter via electron recoil. While it would be interesting to investigate the effect of the fuller sample of DM velocity distributions, this is beyond the scope of the current paper but in any case it is likely to only further strengthen our findings since the fuller sample would likely produce a broadened range of the astrophysical parameters of relevance to this analysis.
In figure 1 we present the local dark matter speed distributions in the Galactic rest frame for the 36 MW-like halos selected from the TNG100 simulations that exhibit MB velocity distribution. As can be seen in the figure, there is quite large halo-to-halo variation in the speed distribution. The SHM peak speeds are in the range of v 0 = [181-238] km/s and the escape speeds are in the range v esc = [435-558] km/s. The peak speeds of the distribution are found by fitting eq. (2.1) while the escape speed is taken to be the highest speed of any particle observed in the torus around the Sun's location. An alternative definition of the latter would be to choose the speed of the particle needed to escape from the potential well of all of the particles in the halo. We have verified that using this definition does not alter the main results of the paper.
The derived ranges are consistent with, though somewhat larger than, those found, for example, in [52]. The relatively large spread in the speed distributions owes to the fact that in this work we find many more halos that fulfill the selection criteria due to using a lower resolution simulation within a larger cosmological volume. The point here is that the IllustrisTNG100 simulations employ a simulation box size of roughly 100 Mpc side length. Previous studies of the dark matter velocity distribution in cosmological simulations either used simulations with higher resolution (e.g. EAGLE simulations with a box of 25 Mpc side length in reference [70]) or used zoom-in simulations where a few halos are preselected and then re-simulated with a higher resolution. Therefore, statistically the larger cosmological volume is assumed the more galaxies it will contain and the more MW-like galaxies can be identified.
Nevertheless, it would be interesting to cross-check our results within the same suite of simulations but using a higher resolution box. This can become possible once the results of TNG50 become publicly available [71].
In order to study the impact of inferred ranges of the speed distributions on the direct DM searches we choose the halos with the minimum and the maximum peak speeds, shown in figure 1 in green and red, respectively. We also define the median values 2 of the peak speed and the escape velocity, the corresponding velocity distribution is shown in figure 1 in blue. We refer to them as MIN, MED and MAX profiles, as they correspond to minimal, median and maximal expected strength of the signal in direct detection experiments. We present the results in table 1. Additionally, for each case we calculate the local dark matter density in the torus region, that we also list in table 1. The resulting local velocity profiles are first of our main results of this work.

JHEP07(2020)081
Name 3 Uncertainties in sub-GeV DM electron recoil As an illustrative and timely example we apply the DM speed distribution profiles found in the previous section to direct detection searches of sub-GeV dark matter. The impact of astrophysical uncertainties on the sub-GeV dark matter range is a largely unexplored area, with most studies so far focusing on WIMP detection through standard nuclear recoil experiments [72][73][74][75][76][77][78]. These studies accommodated a wide range of circular and escape velocities compatible with the calculations given in the Illustris simulation in this work. 3 However, nuclear recoil experiments are insensitive to dark matter in the sub-GeV range.
For such light DM one of the most promising detection techniques is based on electron recoils, which can provide a strong signal through electron ionization in small band-gap semiconductors. For silicon and germanium this gap is typically ∼1 eV, and so allows for detection of dark matter even as light as hundreds of keV. The scattering DM-electron rate in a semiconductor, assuming that DM speed distribution f (v) follows a symmetric MB distribution, 4 is given in events/kg/day by [10] the following: where ρ χ and m χ are the local density and mass of the DM particle, respectively, and m e is the mass of the electron. The crystal form factor, f crystal (q, E e ) contains information about the electronic structure of the material, N cell is the number of unit cells in the crystal target, α is the fine structure constant andσ e parameterizes the strength of the interaction, which in the case of F (q) = 1 is equal to the cross-section for an elastic scattering of DM on a free electron. We also consider the case where F (q) ∝ 1/q 2 which corresponds to electron-DM scattering through a light mediator.

JHEP07(2020)081
To calculate the rate for the detector we must first boost the DM velocity with the average Earth velocity relative to the DM halo, taken to be v E = 240 km/s such that we define g(v) ≡ f (v + v E ). With this consideration, the astrophysical properties of the halo are then encapsulated by the function η (v min ) such that where, through energy conservation, v min is a function of the electron energy E e and momentum transfer q, This determines the minimum speed DM particles required to excite electrons in the crystal, transferring q momentum and energy E e . The stark contrast between electron excitation and nuclear recoil based experiments is the fact that the electron is much lighter and faster than a nucleon. Therefore, the momentum transfer q will be approximately given by the momentum of the electron q ≈ m e v e , where v e is the electron velocity which is typically an order of magnitude larger than the circular velocity of the halo, v 0 . Hence, energy transitions to states above the band gap are sensitive to the distribution around the tail of the DM velocity profile. The interaction between an incoming DM particle and an electron in the crystal is in fact a rather complicated process. It is understood that after the DM particle deposits the initial energy E e to the electron, additional secondary scatterings redistribute this energy between lower-lying electron or hole energy states. Eventually the secondary scatterings fall below a certain energy threshold for electron-hole pair production, ε, and the number of total pairs generated is what is measured. Assuming pair creation is linear in energy, and following the conventions in [10], we define the ionization threshold Q as a function of the electron energy E e , where the energy band gaps E gap for silicon and germanium are 1.11 eV and 0.67 eV, respectively. The energy needed to produce an additional electron-hole pair above the band gap, denoted by ε, is also determined to be 3.6 eV and 2.9 eV for silicon and germanium, respectively. The expression (3.4) takes only integer values for Q since it represents the total number of electron-hole pairs produced for a single DM interaction. This is a key parameter in calculating experimental sensitivity curves for direct detection. For example, the DAMIC experiment that is currently running at SNOLAB [80] expects to reach a threshold as low as 2-3 electrons in a total of 1 kg of detector with its current level of leakage current (produced by thermal excitation in the crystal) demonstrated to be very low in the prototype run [81]. Among other considerations, improvements in leakage current reduction through operation temperature and material quality is imperative in reaching low thresholds. Hence, in this study we examine the consequences of halo DM properties on multiple thresholds.

JHEP07(2020)081
To calculate the direct detection limits for silicon and germanium targets, we employ the publicly available code QEdark [10] in accordance with the plane wave self-consistent field (PWscf) code in Quantum ESPRESSO [82] for the electron wave-functions and energy levels. QEdark is utilized to numerically compute the crystal form factors for both cases, and then subsequently the total scattering rate and sensitivity curves.
To illustrate the impact of different velocity profiles in figure 2 we compare normalized event rates for the MIN (green), MED (blue) and MAX (red) cases, for Ge and Si targets and two different DM mass values. The large difference in the velocity tails of these profiles, parametrized mostly by the value of the v 0 , translates to more extended tails in the electron energy distribution E e . This leads to higher expected rates, with the difference being more pronounced in the higher ionization thresholds. The characteristic shape of the rate histogram for germanium (bottom right panel) is a result of the fastmoving 3d-shell electrons that dominate at around E e 24 eV. Despite interacting with a larger momentum transfer with the DM particle, their contributions are not significant for the smaller thresholds.
The final results for projected 95% C.L. sensitivities for a 1 kg-year exposure with the astrophysical uncertainties taken into account are given in figure 3 for silicon and in figure 4 for germanium targets for F (q) = 1 form factor and in figures 5 and 6 corresponding plots for F (q) ∝ 1/q 2 . Left panels show the projected limits on the cross-sectionσ e and right show the relative size of the uncertainty due to the DM density and velocity profile. In the limit of very large m χ the uncertainty reaches a constant, since the minimal velocity essentially becomes independent of the DM mass, see eq. (3.3). For light DM masses near the detection thresholds the uncertainty in the sensitivity grows significantly where both the circular and escape velocity of the DM profile contribute sizable effects. This shifts the precise position of the minimum reachable mass for a particular pair threshold. The projected sensitivities are calculated under the assumption of no background, with the understanding that it is effectively parameterized by the value of Q that is achievable (for more detailed discussion of the backgrounds see [10]).
What is worth stressing is that: (a) the impact of the velocity profile is quite significant, and is comparable to the change between consecutive electron-hole pair thresholds, and (b) the lower the electron-hole pair threshold is experimentally achievable, the better is not only the overall sensitivity but also the robustness of the detection method.

Discussion and conclusions
In this paper we analyzed the results of the recent IllustrisTNG hydrodynamical simulation and obtained the local DM halo properties. We identified the resulting best fit (MED) as well as -from the point of view of predicted sensitivity of DM direct detection experiments -extreme MIN and MAX profiles. These profiles amount to first of the two main results of this work and can be adopted in updated direct detection analyses. Next we used them to quantify, for the first time, the astrophysical uncertainties of the sub-GeV electron recoil DM searches. Qualitatively, based on experience with nuclear recoil based searches, one would expect that: (a) the impact of the velocity profile is to shift the DM minimal mass threshold, and (b) that the effect is most pronounced for m χ ∼ m e . In many electron recoil based detection techniques, however, because of their sensitivity to the energy deposition in generating the measurable signal even for m χ m e , the velocity profile of DM can play a significant role. We expect this to be the case for not only semiconductor-based detectors, but potentially also detectors using Dirac materials or superconductors. In contrast, experiments measuring only the direct effect of single recoil or absorption, without sensitivity to the exact value of the energy transfer, suffer from much smaller astrophysical uncertainties, more closely resembling the known case of nuclear recoils. Quantitatively, our results show that the uncertainty reaches more than an order of magnitude for two most studied semiconductor targets, silicon and germanium. More  Figure 6. The same as figure 4 but for F (q) ∝ 1/q 2 form factor. precisely, even in the limit of large DM mass (m χ 100 MeV), the sensitivity difference between MIN and MAX velocity profiles reaches factor O(8-16) depending on the target material and more importantly the ionization threshold. This is even for the contactlike interaction form factor F (q) = 1, while other form factors favoring interactions with smaller momentum transfer, lead to larger effect. E.g. for F (q) ∝ 1/q 2 the sensitivity difference is O (12)(13)(14)(15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30). Moreover, for lower DM masses the effect is much larger, due to higher sensitivity of the tail of the velocity distribution.
Therefore, to summarize, the second main result of this work underlines how crucial it is to develop background reduction techniques for experiments based on electron recoils, as not only the sensitivity grows and theoretical uncertainty diminishes, but also the impact of the not yet well known dark matter halo profile is much less severe if one is able to be sensitive to the lowest possible ionization thresholds. Note added. After our paper was submitted to arXiv reference [83] appeared which focused on future sensitivity projections of sub-GeV DD experiments, where the impact of varying the astrophysical parameters of the SHM is also addressed at some level.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.