Dark matter substructures affect dark matter-electron scattering in xenon-based direct detection experiments

Recent sky surveys have discovered a large number of stellar substructures. It is highly likely that there are dark matter (DM) counterparts to these stellar substructures. We examine the implications of DM substructures for electron recoil (ER) direct detection (DD) rates in dual phase xenon experiments. We have utilized the results of the LAMOST survey and considered a few benchmark substructures in our analysis. Assuming that these substructures constitute $\sim 10\%$ of the local DM density, we study the discovery limits of DM-electron scattering cross sections considering one kg-year exposure and 1, 2, and 3 electron thresholds. With this exposure and threshold, it is possible to observe the effect of the considered DM substructure for the currently allowed parameter space. We also explore the sensitivity of these experiments in resolving the DM substructure fraction. For all the considered cases, we observe that DM having mass $\mathcal{O}(10)\,$MeV has a better prospect in resolving substructure fraction as compared to $\mathcal{O}(100)\,$MeV scale DM. We also find that within the currently allowed DM-electron scattering cross-section; these experiments can resolve the substructure fraction (provided it has a non-negligible contribution to the local DM density) with good accuracy for $\mathcal{O}(10)\,$MeV DM mass with one electron threshold.

typically achieved at a lower DM mass. For e.g., assuming a xenon target and momentum independent scattering cross-section, the maximum sensitivity is achieved at ∼ 30 GeV for DM-nuclear scattering and ∼ 200 MeV for DM-electron scattering.
An ambient DM of mass O(10) MeV will have a kinetic energy of the O(10) eV, which is in the ball-park of the atomic ionization energy or the band gap energy of semiconductor. This indicates that a sub-GeV DM can ionize an electron from an atomic shell or facilitate an electron's transition from the valance band to the conduction band. Many experiments like XENON [63], SuperCDMS [64], DarkSide-50 [65,66], DAMIC [67], EDELWEISS [68], SENSEI [69,70], PandaX-II [71] etc. are searching for the signatures of such a phenomenon.
The boundedness of electrons in the target material makes DM-electron scattering events inelastic. The DM velocity required to have a measurable recoil is rather high, which can be found near the tail of the DM velocity distribution (assuming that it has a Maxwell-Boltzmann form). These tails are quite sensitive to the choice of the DM velocity distribution [72][73][74][75]. The present DM velocity distribution depends on the galactic structure formation history. In the well-known paradigm of ΛCDM (Lambda Cold Dark Matter), bottom-up hierarchical structure formation is a generic feature [76][77][78][79][80][81][82]. Larger galaxies are formed from the merger of smaller galaxies (although the merger of similar mass galaxies may also lead to a bigger galaxy [83,84]). The gravitational field of the Milky Way (MW) is non-uniform, and this non-uniformity gives rise to strong tidal forces. When smaller galaxies accrete into the MW galaxy, the gravitational force disrupts these galaxies resulting in tidal stripping of various components (including DM) of these infalling galaxies. For an ancient merger, the DM component will have time to virialize within the MW, which may lead to an isotropic, isothermal DM halo. This scenario is often referred to as the Standard Halo Model (SHM), with the Maxwell-Boltzmann distribution representing the DM distribution. However, for relatively recent mergers, there will not be sufficient time for virialization, resulting in plenty of substructures both in the stellar and in the DM component [85][86][87][88][89][90][91][92][93][94][95][96][97]. The presence of such additional stellar substructures (beyond the MW stars) have been detected by different sky-surveys like Gaia [83,90,[97][98][99][100], SDSS [90], LAMOST [101,102], etc., and have also been predicted in various N-body simulations [103][104][105][106][107][108][109][110][111].
Since these stellar substructures arise from merged galaxies, a DM counterpart must be associated with them too (because the DM is also present in the accreted galaxies before their merger). Whether DM would follow stellar distribution or not is a matter of debate. For example, the celestial part of the Sagittarius stream might not substantially overlap with the Solar neighborhood. However, the extended DM counterpart may overlap with our local position [112]. The similarities between DM and stellar distributions in debris flow have been pointed out in Refs. [89,113]. The dwarf spheroidals, which give rise to the S2-stream, are believed to have similar DM and stellar shape [114] before they merged with MW. Therefore the resemblance between stellar and DM substructures is not settled yet; more dedicated studies are needed to understand this. However, the presence of this DM might manifest in the local DM density and velocity distribution: this will result in a difference of the velocity distribution from the normal MB distribution with cut off at the galactic escape velocity [115,116]. DM DD rate is strongly dependent on the local velocity distribution of DM , and a different DM velocity distribution can result in a large change in our theoretical expectations. The effects of these substructures have been extensively studied in the literature in the context of typical NR DD experiments [114,[147][148][149][150][151][152][153][154][155][156][157][158][159][160]. This paper aims to study the effect of these DM substructures in the ER DM DD experiments assuming xenon-based detectors. Such a study has been conducted for semiconductor target material in Ref. [73]. It was shown in Ref. [75] that the effect of such astrophysical uncertainties is quite prominent for xenon targets. Further, in large regions of the DM parameter space, the sensitivity of xenon targets is a few orders of magnitude stronger than those from semiconductor-based experiments [63,[69][70][71] implying that xenon detectors will probably play a big role in discovering DM-electron scattering. These facts motivate our detailed study in this manuscript, where we highlight the importance of considering DM substructures while searching for DM-electron scattering.
It has been argued in Refs. [83, 90, 92-94, 98, 99, 161] that there are plenty of stellar substructures in the local halo. We utilize the results of the LAMOST survey [101] to present the effect of the DM substructure [94] in DM ER experiments. Without a loss of generality, we demonstrate our results by choosing a few benchmark substructures. We expect broadly similar results for other relevant substructures. In addition, our formalism will be useful for future analysis of DM ER experiments for xenon-based targets. Currently we do not understand how much of these substructures contributes to the local DM density. We have chosen a few benchmark values of DM substructure contributions to the local DM density, namely 100%, 20%, and 10% and presented our results. Our choices are motivated by Ref. [96] which states that stellar substructures near the Sun may constitute 20% of the stellar halo. We also consider the forecast of xenon targets in resolving the fraction of DM substructures components for a few benchmark choices of the DM parameter space.
The rest of the paper is organized as follows. In Sec. 2, we briefly review the DMelectron scattering in xenon-based detectors. In Sec. 3, we describe DM substructures that we have considered in our analysis. In Sec. 4, we present our results along with the statistical methodology, and conclude in Sec. 5.

DM-electron scattering at xenon
If the ambient DM particle scatters off an electron of xenon, DM may transfer its kinetic energy to the electrons, leading to free electrons. For example, a non-relativistically moving ambient DM of mass ∼ 100 MeV will have kinetic energy ∼ 50 eV (in the Solar system), which is in the ballpark of the electron ionization energy of xenon.
In a two-phase xenon time projection chamber, DM particles interact with the liquid Xe target material, and depending on interaction type (electronic or nuclear), the signal topologies are different. For DM-nuclear interaction, the deposited DM energy produces excited atoms, electron-ion pairs, and some non-observable heat. Some free electrons recombine with ionized atoms to generate more excited atoms. Essentially both the direct and excited states produced by electron-ion recombination make a characteristic scintillation light. This prompt scintillation light, known as S1, is detected in photomultiplier tubes (PMTs) immersed in the liquid Xe at the bottom. Due to an external electric field, the remaining electrons drift through liquid xenon and cross the liquid and gaseous interface, producing proportional scintillation in the upper PMTs. This signal is known as S2. For the ER interactions, almost all the ionized electrons are collected at the upper PMTs through scintillation, producing a dominant S2 signal with a subdominant S1 signal. Hence ER interactions manifest through a large S2/S1 ratio compared to the NR case [162].
Let us consider a DM particle of mass m χ and velocity v scattering off an electron in the xenon atom. Energy conservation implies [60] v where v min is the minimum DM velocity required to get an ER of ∆E e , and q is the momentum transfer to the electron. Note that ∆E e must be greater than the ionization energy of the corresponding shell E n,l to have an observable recoil E e , i.e., ∆E e = E n,l +E e . The differential DM-electron scattering event rate can be written as [57] where N T is the number of electrons in the target, ρ χ denotes the local DM density, and DM-electron reduced mass is represented by µ χe . DM-electron scattering cross section for a reference momentum transfer, namely q = αm e , is indicated byσ e . The DM form factor, F DM (q), takes care of the momentum dependency in the cross-section. The ionization form factor is represented by f n,l ion with n and l being the principal and angular momentum quantum number, respectively. The recoil momentum is denoted by k = √ 2m e E e . The time dependency of the recoil signal is described through t. The quantity η, also called the mean inverse speed, depends on the i th DM velocity distribution as where f i lab is the DM velocity distribution at the detector's rest frame in the location of the Earth for the i th DM component (which contributes to the DM velocity distribution). The latter can be obtained by boosting the galactic rest frame DM velocity distribution where v E is the Earth's velocity in the galactic rest frame: Here v LSR is the velocity of the local standard of rest (LSR), v pec is the peculiar velocity of the Sun with respect to the LSR. Conventionally these are expressed in galactic rectangular co-ordinate and expressed as v LSR = (0, v 0 , 0), v pec = (11.1 ± 1.5, 12.2 ± 2, 7.3 ± 1) km/s [163]. Following Refs. [75,156], throughout the paper we fix v 0 = 233 km/s. The uncertainties associated with v 0 and other astrophysical parameters have been studied in Refs. [72,74,75] in the context of ER (see Ref. [164] for halo independent analysis). The time-dependent Earth's velocity is represented by u E (t) which leads to the well-known annual modulation of the signal. The expression for u E (t) can be found in [165]. The differential event rate given in Eq. (2.2) can be divided into three parts. The particle physics input is indicated byσ e and F DM . Throughout our analysis, we will do a model-independent analysis with two choices of F DM : 1 and 1/q 2 , which appears in large classes of particle physics model [166][167][168][169][170][171][172][173]. We will present the results of F DM = 1 in the main text and that of F DM = 1/q 2 in the appendix. The atomic physics part symbolized by f n,l ion signify ionization probability. The numerical values of the f n,l ion is adopted from QEdark [55,57,174]. The local DM density and η constitute the astrophysical inputs.
The galactic DM velocity distribution is traditionally assumed to be a Maxwell-Boltzmann (MB) distribution truncated at the galactic escape velocity (v esc ) and erf is the error function. Throughout the discussion the galactic escape velocity (v esc ) has been fixed to 528 km/s [156,175]. While the MB distribution may describe the DM velocity distribution which is in equilibrium (hydrodynamical simulations indicate that MB distributions may not adequately describe the velocity distribution of the smooth DM halo component), the equilibration condition will not be met for relatively recent mergers of the MW with other galaxies. These recent mergers will have unique signatures, both in velocity and position space, called substructures. The existence of these substructures is also observed in various N-body simulations. When a galaxy accretes into the Milky Way, the stellar component of the accreted galaxy carries several tell-tale signatures: stellar streams, stellar shards, and stellar debris flow [85][86][87][88][89][90][91][92][93][94][95][96][97].
The recent results of various surveys like Gaia, SDSS, and LAMOST indeed indicate the presence of these stellar substructures. Combining the effect of the substructure with the SHM, we get total average inverse speed as where f ζ i lab (v, t) refers to the substructure velocity distribution (discussed in Sec. 3) and δ represents the fractional contribution that the corresponding component constitutes to the local density of DM. 2 In what follows, we will consider the effect of these substructures in DM velocity distribution and the ER DD rate in liquid xenon experiments.

DM substructures
This section discusses the benchmark DM substructures that we have studied in this work. We have utilized the results of Ref. [94] where the stellar substructure is obtained using  [94]. The DTG from which substructures are identified has also been specified.
the star catalog of LAMOST DR3 [101]. We choose a few representative substructures to present our results. For clarity, we also mention the name of the associated dynamically tagged groups (DTG) with the relevant substructures [94]. The details of these substructures are summarised in Table 1. We emphasize that the chosen substructures are for illustrative purposes only. Further research is required in order to understand the DM content of various substructures and whether the substructure DM profile coincides with the Solar circle. Whether the corresponding DM substructure will follow the same velocity distribution as the stellar substructure or not is currently not understood. Using Via Lactea II high-resolution N -body simulation, it has been shown that DM debris flows closely follows their stellar counterpart [89,113]. However, the same is not valid for Sagittarius stream [112]. Nevertheless, we will assume that the velocity distributions of the substructures follow that of the corresponding stellar components. This assumption can be confirmed or refuted by future research. However, the broad conclusion (like the change in the event rate and subsequently in the discovery limit due to DM substructures) of this study will hold. We note that the substructures we have considered in this paper have similarities with previous considerations [73,114]. For instance, the Helmi substructure is analogous to S2substructure [109]. The velocity properties of the Nyx substructure are somewhat similar to the prograde (Pg) stream and are expected to arise from the same Splashed Disk event [94]. 3 Some of the considered substructures are also found in Gaia DR3 data at the Solar neighborhood [95,96,161].
The mean stellar velocities and the diagonal values of the stellar velocity dispersions are given in Table 1. In general, DM substructures will have a different velocity distribution than the virialized component (SHM), which will dramatically impact the ER distribution. The galactic velocity distribution for each of the substructures (referred to by ζ i ) can be written as [73,114] Figure 1: We show the lab frame DM speed distributions for all the considered astrophysical components. The red, brown, orange, olive, blue, purple, and pink solid lines represent the HelmiDTG1, HelmiDTG3, PolarDTG11, PgDTG2, Sausage, RgDTG28, and Sequoia, respectively. We also display the same for SHM by the solid black line. The vertical dashed line shows required v min for m χ = 100 MeV with E e = 20 eV, q = 25 keV, and 5p 6 shell.
where σ ζ i is the velocity dispersion matrix, assumed to be diagonal with the values given in Table 1 and det σ ζ i is the determinant of the dispersion matrix. The mean velocities of the substructures in the galactic frame are expressed by µ ζ i which are non-zero in contrast to the SHM case, as indicated in Table 1. The normalization constant N ζ i esc is calculated numerically. The step function represents the cut-off at the galactic escape velocity, although the substructures' velocity distributions are likely to peak at smaller velocities. Therefore this cut-off will have a numerically insignificant effect. The index ζ i refers only to the substructure, whereas i includes both the substructures and SHM.
Assuming Eq. (3.1) as the galactic velocity distributions for the DM substructures, we display the corresponding lab frame speed distributions, Fig. 1. Except for the modulation signature (discussed in Sec. 4.3), we fix the Earth's velocity to v E = (39.7, 243.2, 16.4) km/s to economize the computation. This value of v E is attained during first week of March when the Earth's velocity is roughly equal to its average velocity. For Sequoia we have explicitly checked that taking the exact yearly average rate would lead to less than ∼ 6% change in the discovery limit. Given the poor knowledge of DM substructure fraction, we ignored this 10% effect. The general trend we observe is that the substructures which peak at larger values of v have negative µ φ . Since the Earth moves with high positive rotational velocity ∼ 250 km/s, substructures with negative µ φ , will hit the Solar system with larger velocities. On the other hand, substructures having large positive µ φ co-rotate with the Earth, leading to f i lab (v) peaking at smaller velocities. This has been displayed in Fig. 1, where the Helmi streams having larger values of µ φ peak at relatively smaller velocities, whereas Sequoia having a negative µ φ peaks at the higher velocity. We also display the velocity distribution of SHM by the solid black line. For reference we show the required v min = 428.7 km/s to obtain a recoil of 20 eV with momentum transfer 25 keV and 5p 6 shell for DM mass 100 MeV by the vertical black dashed line.
Given these velocity distributions, we turn to the discussion of the mean inverse speed η i (v min ) (using Eq. (2.3)) of each of the astrophysical components. The values of η i (v min ) as a function of v min are depicted in Fig. 2. The vertical black dashed line indicates values of η i (v min ) for v min = 428.7 km/s. Expectedly, η i (v min ) are monotonically decreasing function of v min , which can be understood from the integration over velocity starting from v min . The maximum values of η i (v min ), i.e., η i (0) is larger for the distributions which peak at lower velocities because the mean inverse speed is inversely proportional to the most probable speed (the speed at which velocity distribution attains maximum value) of the distribution. Hence in Fig. 2, we observe maximum and minimum η i (0) for HelmiDTG3 and Sequoia respectively. For the other distributions, η i (0) lie within the same of HelmiDTG3 and Sequoia. The flatness of η i (v min ) for Sequoia up to a large value of v min as compared to other distributions is also a manifestation of the higher most probable speed of Sequoia. This indicates the extent to which v min is supported by the distribution. It should also be noted that the flatness of η i (v min ) is also sensitive to the choice of the velocity dispersion.

DM-electron scattering at xenon: effect of substructure
In this section, we discuss the effect of the substructures on the DM-electron scattering rate for liquid xenon experiments. For F DM (q) = 1, the constraint on the DM-electron scattering cross-section from the xenon detectors dominate when DM mass is 50 MeV. Xenon experiments may have a better prospect of discovering DM-electron scattering, and it is essential that we study this prospect thoroughly. Our work outlines the theory effort toward answering this important question.
Following Ref. [57], we convert the ER energy (E e ) to number of electrons (n e ). DMelectron scattering would produce n e number of observable electrons, unobservable photons, and heat. Some primary electrons would recombine with secondary ions with probability f R . Further, each recoiling electron of energy E e will give rise to additional secondary n (1) e = Floor[E e /W ] quanta (photon or electron). The average energy required to create a single quanta is W . Moreover, the scattering process can also lead to the ionization of electrons from the inner shell, which would de-excite by releasing a photon. These photons may also create secondary quanta, n 2) which will give the differential event rate as a function of number of produced electrons. Our paper does not consider uncertainties associated with W, f e , and f R .
In Fig. 3, we show the differential event rate as a function of n e for m χ = 100 MeV, σ e = 10 −40 cm 2 , and 1 kg-year exposure. For each event rate, we have assumed that the corresponding astrophysical component (SHM or substructures) constitutes 100% of the local DM density. For m χ = 100 MeV with typical momentum transfer of O(10) keV, to obtain a measurable recoil the required minimum DM velocity should be around 500 km/s. Hence, the tail of η i (v min ) dominantly contributes to the recoil rate. Evidently the substructures having the largest value of η i (v min ) near v min ∼ 500 km/s give rise to a larger event rate.

Neutrino background
The scattering of neutrinos with electron/ nucleon may also give rise to ionization signals in low-threshold DD experiments. Other background sources like radioactive background, Cherenkov radiation, etc. which can potentially mimic a DM signal [177]. The experimental collaborations confront and beat these non-neutrino backgrounds using various experimental techniques to isolate a potential DM signal. However, the neutrinos are an irreducible background that can not be removed by using shielding, purified detector material, and other experimental techniques. Because of this, we have taken neutrinos as the only source of background in our analysis. If other non-neutrino backgrounds are found in the data-set, then our results will degrade proportionally.
It has been argued in Refs. [178,179] that Solar neutrinos are the main source of back- ground for sub-GeV DM-electron scattering. 4 Neutrino-electron elastic scattering is the dominant contribution of background events for rather large recoil energy (∼ 10 5 eV). Instead, coherent neutrino-nucleon scattering may produce small ionization, which would be the dominant source of background in our consideration. The neutrino-nucleon scattering event rate is [178,181] where N T , M , and T are the number of target nuclei per unit mass, total mass, and time respectively. The minimum neutrino energy to produce a nuclear recoil of energy E NR is expressed by E min ν = m N E R /2. The differential coherent neutrino nucleon cross section and the differential neutrino flux are denoted by dσ/dE NR and dφ ν /dE ν respectively [178,182]. We have utilized low, fiducial, and high ionization models given in Ref. [178] to obtain number of electron n e for a particular nuclear recoil energy. The corresponding neutrino-induced event rate for fiducial model is displayed in Fig. 3 by the grey dashed lines. 5 The grey shaded regions represent variation in the event rate for high and low ionization models of n e [178]. Since there is a difference between three ionization models in the low n e /energy bins, hence we observe a large change in the differential event rates at those bins. The discovery limits for low and high ionization models is given in appendix C. 4 See Refs. [178,180] or discussion related to the prospect of these detectors in probing beyond SM interactions of neutrino. 5 We note that there is a factor ∼ 3 difference in the event rate between our result and Ref. [178].
For one electron threshold, the impact of the ionization model uncertainty leads to less than a factor of 3 change in the discovery limits.

Statistical methodology
In this section, we discuss the statistical procedure to obtain the discovery limit for DMelectron scattering in the presence of substructures for liquid xenon experiments. We have employed the profile likelihood ratio test [183] withσ e and substructure fraction (δ) as the signal parameters of interest. In the following, we briefly discuss this procedure. The binned likelihood for the background and signal model (M), is given by Here θ = (m χ ,σ e , δ, Φ) and D(θ) is the Asimov data set. The number of energy bin is represented by n bins . The Poisson probability (P) at the i-th bin is calculated using observed N i obs and the expected number of events. The expected number of events is the addition of DM events (N i χ ) and the sum of neutrino events (n i ν ) for all the neutrino components (n ν ). The Gaussian function (G(Φ j )) takes care of the uncertainty in the neutrino fluxes (Φ j ) with mean values and standard deviation given in [178,182].
Depending on the choice of the analysis, we vary one of the signal parameters (eitherσ e or δ), treating the other one as a nuisance parameter. We treatσ e as the signal parameter for the discovery reach. Therefore, the profile likelihood ratio test statistic, which compares the background-only hypothesis (M 0 ) with the background and signal model (M), is given by [73,183,184] where λ contains the nuisance parameters, i.e., δ and Φ j in this case. Using Wilks' theorem, it can be shown that the ratio in Eq. (4.3) follows a χ 2 distribution with one degree of freedom [73,183]. Thus, the significance of rejecting the background-only hypothesis is given by √ q 0 -σ. In this paper, we present all the discovery limits at the 90% confidence level (CL). We obtain the discovery limits utilizing Asimov data set which assumes that the number of observed events is same as the expected events. However in a real experimental data-set, this will not be true and in that case one should treat experimental data as the observed events. Then it would be possible to constrainσ e assuming a value of substructure fraction. We consider δ as the signal parameter for the prospective detection of DM substructure fraction. The corresponding profile likelihood ratio test to distinguish two neighboring points δ 1 and δ 2 can be written as [73] This profile likelihood ratio is employed to reject the null hypothesis, which is that two neighboring points δ 1 and δ 2 are indistinguishable at 68% CL. Both for Eqns. (4.3) and  Figure 4: Discovery limits at the 90% CL for different DM substructures considered in this paper, assuming that each astrophysical component accounts for 100% of the local DM density. The discovery limits are obtained using one kg-year exposure, one electron threshold, and F DM = 1.
(4.4) we utilized Asimov data set [183] to obtain the likelihood ratio test. In this scenario, artificial data is generated using the model's parameters (in our case M). Then the expectation is that the number of observed events (N obs ) should be equal to the number of the expected event (N exp ). For a sufficiently large number of observations, the value of the profile likelihood ratio test approaches the median value. Compared to the Monte Carlo simulation, the Asimov data set scenario is computationally more economical while acquiring accurate results. For the 68% and 90% CL limit the required q 0 's are 0.99 and 2.71 respectively. For a fixed m χ and δ, the 90% CL discovery limit is obtained by changinḡ σ e in Eq. (4.3) until the required q 0 (2.71) is achieved. The 68% CL contours in resolving substructure fraction are estimated using Eq. (4.4). In this case for a fixed values of m χ , σ e , and δ 1 , we iterate over δ 2 until the required q 0 (= 0.99) is attained.

Results
Here we will present the results using the statistical analysis discussed in the previous subsection. The three parameters of interest are DM mass (m χ ), DM-electron cross section (σ e ), and the DM substructure fraction (δ). Given that DM has to be massive, we present our results through two possible choices, keeping one of the other two parameters to a fixed value. In the first part, the results are presented through the discovery limit, which is depicted in DM mass and DM-electron cross-section plane keeping a fixed DM substructure fraction. In the other case, considering a fixed DM-electron cross-section, we present the forecast of the xenon experiments to resolve the substructure fraction for a few benchmark choices of DM particle masses. In Fig. 4, we present the sensitivity to DM-electron cross-sections for each of the substructures considered in this paper, assuming that the corresponding substructure constitutes 100% of the local DM density. In Fig. 4, each line represents the minimum DM-electron cross section required to observe the effect of the corresponding substructure in a liquid xenon detector with 1 kg-year exposure and one electron threshold. The discovery limits for two and three electron thresholds are given in appendix B. The different discovery limits for different substructures are the implication of non-identical most probable speed. The tail of the DM velocity distribution will be more populous for the substructure having a relatively larger most probable speed. Therefore a sizable number of DM particles will be available to interact with the target electrons. This leads to a larger event rate, as has been depicted in Fig. 3, where for a fixed DM-electron cross section among the considered DM substructures, we obtain the minimum and the maximum number of events for HelmiDTG3 (lowest most probable speed, see Fig. 1a) and Sequoia (highest most probable speed, see Fig. 1b) respectively. Owing to this, the DMelectron cross-section that can be probed for HelmiDTG3 is the largest, whereas the same for Sequoia is the lowest. The event rates and subsequently the discovery limits lie between HelmiDTG3 and Sequoia for the other considered substructures. The light grey shaded region demonstrates the constraint from the ionization signals in the XENON1T experiment [63], which is the most stringent current DD constraint for the parameter space shown in the plot. For reference, we have also shown the discovery limit for the SHM with the solid black line.
In reality, these substructures may not contribute 100% to the local DM density. Therefore, we choose two benchmark values of δ, namely δ = 0.1 and δ = 0.2 (shown in Fig. 5).  Figure 6: 68% CL contours in m χ − δ plane, representing forecast in resolving DM substructure fraction for a few benchmark points. We have assumed 1 kg-year exposure, one electron threshold,σ e = 10 −40 cm 2 , and F DM = 1. In the left and right panels, we show the results for HelmiDTG3 and Sequoia, respectively.
Further, as mentioned above, we have only considered two substructures, HelmiDTG3 and Sequoia, which lie at two extreme ends. SHM constitutes the rest of the local DM density for the combined DM distribution. If the discovery limit for a particular substructure (with δ = 1) is larger compared to SHM, then the same for the combined DM distribution will lie above the SHM limit. This effect would be more pronounced upon increasing δ. In Fig. 5, the combined discovery limit for HelmiDTG3 and Sequoia is displayed by brown and purple lines, respectively. Notably, brown and purple lines lie above and below the SHM scenario. Upon increasing the δ, we observe more deviation from SHM. Importantly, it is still possible to see the effect of these substructures in liquid xenon experiments with this kind of realistic choice of δ. Next, we turn into the discussion of resolving substructure fractions in liquid xenon experiments. Again, we have restricted ourselves to HelmiDTG3 and Sequoia among the considered substructures as these two reside in the extreme ends. The sensitivity in resolving DM substructure at 68% CL is displayed in Fig. 6 for 1 kgyear exposure, one electron threshold, andσ e = 10 −40 cm 2 , with a few benchmark points. Generically, we observe a better resolution for low DM mass. Comparing Figs. 6a and 6b one can see that we will determine the substructure fraction more accurately for Sequoia compared to HelmiDTG3. This is due to Sequoia's large most probable velocity, which leads to a substantial number of DM-electron scattering events. Generically, it is possible to measure the substructure fraction more accurately, which is moving with a higher most probable speed. For δ = 0.1, with the considered exposure, threshold, andσ e , it is difficult to conclude whether the substructure is contributing to local DM density. Interestingly, for DM mass ∼ 50 MeV, and δ = 0.4, xenon target electron scattering experiments can resolve the substructure fraction with ∼ 50% accuracy. Moreover, the structures of the contours can be understood from Eq. (2.1) and from Fig. 2. Both for the lower and higher DM masses, the inclination of the contours is reversed as we compare HelmiDTG3 with Sequoia. For low DM masses, v min is larger (for fixed q and E e from Eq. (2.1)), therefore it is the tail of the distribution which is contributing to η i . Thus for HelmiDTG3 with low mass DM, an increment in δ will reduce the combined value of η. This reduction could be compensated by increasing DM mass for a fixed number of observed events. This results in a slightly tilted contour towards a higher DM mass. Whereas for higher DM mass (thus smaller v min ), the maximum value of η determines the orientation of contours. For Sequoia, the maximum value of η i is less than that of HelmiDTG3. Hence increasing δ for the former will reduce combined η, which can be elevated by reducing v min , i.e., by increasing DM mass.
We have not discussed a distinctive feature of the DM DD signal: annual modulation [54], where the signal event rates vary with the time of the year in a specified manner. Due to the rotation of the Sun around the MW, there will be DM wind in the Solar rest frame. Due to the Earth's rotational motion around the Sun, the event rate will vary with time. For the SHM, the event rate will be larger (smaller) when the Sun and the Earth travel in opposite (same) directions, respectively. Due to this distinctive feature, which the background cannot mimic, annual modulation events are expected to be less dependent on the background reductions and identifications.
Unlike non-modulation case here we take into account variation of v E over time. The main task for the modulation discovery limit would be to evaluate the event rate against both time and energy (N tim ). For a particular energy bin, we obtain modulation events (N mod ), by subtracting each time bin events from average time bin events (N avg ). We do the same exercise for all the energy bins. The corresponding likelihood can be obtained by taking the difference of their individual Poisson distributions, referred to as the Skellam distribution [185] L(D(θ)|M(θ)) ≡ nt j=1 n bins i=1 e −(N tim (i,j)+Navg(i)) N tim (i, j) N avg (i) where, i and j represent each time and energy bin and n t is the total number of time bins. The modified Bessel function of the first kind is denoted by I N mod . We utilized Eq. (4.5) to obtain test statistics (given in Eq. (4.3)) and subsequently the discovery limit. Following this prescription [73], we find that the modulation discovery limit is weaker than the non-modulation counterpart. For example, with SHM or Sequoia, we observed that the modulation discovery reaches are weaker by a factor ∼ 10 − 100.

Conclusions
The presence of DM in the Universe is well established. Many attempts have been made to discover the connection between DM and SM states. Among them, DD experiments look for the scattering signatures of DM and visible states. There has been a growing interest in the search for light DM (masses 1 GeV) through DD. Ambient non-relativistic DM having mass in the sub-GeV range can not impart sufficient energy to produce a measurable recoil in the typical nuclear recoil DD experiments. Electron, being a light particle, can be an excellent target in detecting such light DM. Many target materials have been considered to identify electronic excitation by the scattering of ambient DM. DM velocity distribution is an integral part of calculating the event rate or the exclusion limit of the DD experiments. DM is also an intrinsic part of structure formation; the history of galaxy formation influences its velocity distribution. While it is difficult to track the velocity distribution of DM, however, it may be manifested through stellar distribution. Surveys like Gaia, SDSS, LAMOST, etc., have made unprecedented progress mapping these stellar distributions. These data reveal the presence of stellar clumps and substructures. It is highly likely that there is a DM counterpart to these stellar substructures, called DM substructure. This paper investigates the prospects of detecting these substructures in low threshold DM DD experiments through elastic DM-electron scattering. Specifically, we have explored the prospect of xenon targets experiments in deciphering this. Note that compared to semiconductor targets experiments (like SENSEI), the xenon targets experiments have better sensitivity in the DM mass range of O(100) MeV.
We utilize the results of the LAMOST survey and choose a few benchmark DM substructures. We emphasize that there is no definite proof of the existence of the DM counterpart to the detected stellar substructures. However, it is likely that they exist. If these DM substructures overlap with the Earth's position, then we can observe the imprint of the same in xenon targets experiments through DM-electron scattering. We find that if the substructure constitutes 10% of the local DM density, then there is a possibility to observe the effect of the substructures in xenon target experiments with the currently allowed DM particle properties. We have also explored the forecast of xenon experiments in resolving the DM substructure fraction. We find that the uncertainty in resolving DM substructure fraction is considerable for higher DM mass compared to lower DM mass. For example, with m χ = 50 MeV,σ e = 10 −40 cm 2 , and one electron threshold in xenon experiments, we can resolve the substructure fraction to ∼ 50% accuracy provided δ ∼ 0.4. The discovery limit and resolving DM substructure fraction are mainly regulated by the most probable velocity of the corresponding velocity distribution. Given this correlation between DD rates and DM velocity distributions, a more detailed understanding of DM substructure is required. High-resolution cosmological simulations and near-future observations will play a crucial role in understanding this. We encourage the experimentalists to continue their excellent work in improving their detector sensitivity so that we are sensitive to such a signal. Our work shows that by pursuing this technique, we will be able to know more about the particle physics and astrophysics of DM and maybe even discover it.

A Event rate
In this appendix we provide event rate for DM mass 50 MeV. This is displayed in Fig. 7.

B Discovery limits for two and three electron threshold
Throughout the main text, we have considered the reach of the xenon experiments for one kg-year exposure and one electron threshold with F DM = 1. Here we present the discovery limit with two and three electron thresholds for δ = 0.1. The results are depicted in Figs. 8a and 8b. With higher thresholds, the expected event numbers decrease; thus, the required cross-section to see the possible effect of the substructure increases. Further, the lowest possible DM mass that can be probed also increases.

C Variation in the discovery limits
As discussed in Sec. 4.1, background event rate from neutrino may change depending on the ionization model. In this appendix, we present the discovery limit for high and low ionization efficiencies models for n e [178]. We display the result in Fig. 9. For each of the substructures, solid lines represent discovery limits for fiducial ionization model and shaded bands show the corresponding uncertainties associated with the ionization models.   Figure 9: Variation in the discovery limits due to high and low ionization efficiency models for n e [178] assuming one electron threshold. For each of the substructures the corresponding light shaded bands represent the uncertainties that may arise from different ionization models of n e . Each of the substructures contribute to 100% to local DM density.

D Momentum dependent DM-electron scattering
In this appendix, we present the discovery limits of the momentum-dependent DM-electron scattering, namely F DM = α 2 m 2 e /q 2 for the considered DM substructures. In this case, we also observe a similar tendency, except that the minimum required DM-electron cross  Figure 10: Discovery limits at the 90% CL for the considered DM distribution with F DM = α 2 m 2 e /q 2 . We have assumed one kg-year exposure and one electron threshold to obtain the discovery limit. The Xenon10 [57] and SENSEI@MINOS [70] limit has also been shown by the grey and blue shaded regions.
section for the discovery of the substructures is larger than the same of F DM = 1. This is displayed in Fig. 10.