Implications of the Gaia Sausage for Dark Matter Nuclear Interactions

The advent of the Gaia era has led to potentially revolutionary understanding of dark matter (DM) dynamics in our galaxy, which has important consequences for direct detection (DD) experiments. In this paper, we study how the empirical DM velocity distribution inferred from Gaia-Sausage, a dominant substructure in the solar neighborhood, affects the interpretation of DD data. We survey different classes of operators in the non-relativistic effective field theory that could arise from several relativistic benchmark models and emphasize that the Gaia velocity distribution could modify both the total number of events as well as the shape of the differential recoil spectra, the two primary observables in DD experiments. Employing the euclideanized signal method, we investigate the effects of the Gaia distribution on reconstructing DM model parameters and identifying the correct DM model given a positive signal at future DD experiments. We find that for light DM with mass ~10 GeV, the Gaia distribution poses an additional challenge for characterizing DM interactions with ordinary matter, which may be addressed by combining complementary DD experiments with different targets and lowering the detection threshold. Meanwhile, for heavy DM with mass above ~30 GeV, depending on the type of DM model, there could be a (moderate) improvement in the sensitivity at future DD experiments.

The advent of the Gaia era has led to potentially revolutionary understanding of dark matter (DM) dynamics in our galaxy, which has important consequences for direct detection (DD) experiments. In this paper, we study how the empirical DM velocity distribution inferred from Gaia-Sausage, a dominant substructure in the solar neighborhood, affects the interpretation of DD data. We survey different classes of operators in the non-relativistic effective field theory that could arise from several relativistic benchmark models and emphasize that the Gaia velocity distribution could modify both the total number of events as well as the shape of the differential recoil spectra, the two primary observables in DD experiments. Employing the euclideanized signal method, we investigate the effects of the Gaia distribution on reconstructing DM model parameters and identifying the correct DM model given a positive signal at future DD experiments. We find that for light DM with mass ∼ 10 GeV, the Gaia distribution poses an additional challenge for characterizing DM interactions with ordinary matter, which may be addressed by combining complementary DD experiments with different targets and lowering the detection threshold. Meanwhile, for heavy DM with mass above ∼ 30 GeV, depending on the type of DM model, there could be a (moderate) improvement in the sensitivity at future DD experiments.

I. INTRODUCTION
Confirming the existence of dark matter (DM) through a variety of cosmological and astrophysical observations has been one of the major successes of 20 th century physics. At the same time, questions regarding the particle nature of DM and its interactions with ordinary matter beyond gravity remain unresolved. Fortunately, there is a vibrant research program that seeks to answer these questions on the experimental and observational frontiers. A leading probe in the hunt for DM is direct detection (DD) experiments, which look for signals from DM particles scattering in underground detectors. Although there have been no statistically significant detections of non-background events so far, next-generation experiments such as Lux-Zeplin (LZ) [1], Xenon-nT [2], PandaX-xT [3], SuperCDMS SNOLAB [4], Damic-M [5], Darwin [6] and Darkside-20k [7] serve as promising avenues not just for DM discovery, but, as we argue below, for reconstructing its astroparticle properties as well.
The main physical observable in a DD experiment is the differential recoil spectrum, typically quoted as a function of the primary scintillation signal. Interestingly, modeling the DM recoil spectra at DD experiments relies on several independent aspects of its phenomenology. More specifically, DD experiments probe a combination of three important DM properties: its mass, interaction type (with nucleus and/or electrons) which we will refer to as the model, and its astrophysical distribution, namely through its density and velocity distribution in the solar neighborhood. From the perspective of statistical inference, this results in a three-fold degeneracy depicted in Fig. 1. As corollary, determining the local astrophysical properties of DM precisely will be crucial in reconstructing its particle physics properties.
The Standard Halo Model (SHM) of DM velocity distribution [8] has been the cornerstone of DD analyses since it was proposed nearly three decades ago. The SHM follows from modeling the Milky Way (MW) as an isotropic, isothermal halo in equilibrium formed through the virialization of multiple subhalo merger residue. However, over the years, observations [9][10][11][12] and results from N -body simulations [13][14][15][16][17] have suggested the presence of diverse stellar and DM substructures from recent mergers, challenging the MW's steady state characterization. There have also been attempts to semi-analytically model the local DM velocity distribution using form of kinematic data [18][19][20][21][22][23][24], but these rely on additional, and potentially restrictive, assumptions about the structure of the MW.
On the other hand, astrometric data released by the European Space Agency's Gaia mission [25,26] presents a unique opportunity to study the MW's accretion history, and take first steps toward an empirical determination of the DM phase-space distribution. Even with a subset of the full data, a few groups [27][28][29][30] have reported evidence of tidal debris from a dominant merger in the solar neighborhood, the so-called Gaia Sausage or Gaia Enceladus, with very different kinematics compared to the old, virialized stellar population. While the full implication of stellar data for the astrophysical properties of DM will take decades to analyse, pioneering work by ref. [31] used a mixture model analysis to characterize the local DM velocity distribution. They used the kinematics of MW halo stars in the cross-matched Sloan Digital Sky Survey (SDSS)-Gaia data set as tracers for the DM velocity, and validated their analysis [32] with the Fire-2 cosmological zoom-in simulation. There has also been growing interest in studying the effect of other phase-space substructures in the solar neighborhood discovered using Gaia data, in particular for the retrograde S1 [33], and the prograde Nyx [34,35] streams. Stellar streams appear as a coherently moving group of stars resulting from the tidal debris of a galaxy localized in both position and velocity space. An important consideration while modeling the DM signal from streams is its non-zero mean velocity in the galactocentric (GC) frame. Thus, presence of a significant DM fraction in stream(s) could result in a very different annual modulation signature compared to the SHM. With this motivation, ref. [36] revisited the DM interpretation of the latest DAMA data with the S1 stream, and found that absence of a DM signal at other experiments rules out the preferred DM parameter space of DAMA even if 100% of the local DM was present in such a stream. Another promising avenue to look for interesting signatures of DM substructure are axion searches and directional detection experiments as discussed in refs. [37][38][39]. The main drawback of the aforementioned analyses is the underlying assumption of a near perfect stellar-DM velocity correlation. In fact, as illustrated in ref. [32] (see top panel of Fig. 7 for example), stellar streams turn out to be poor tracers of the DM velocity, since the tidal debris in the stream hasn't had enough time to completely mix with the halo. Thus we only focus on the DM velocity distribution associated with the Gaia Sausage in our work, and ignore streams until the correlation between their stellar and DM phase space distributions is better established.
In this work, we perform a systematic study of how the DM velocity distribution affects the reconstruction of DM model parameters at current or near-future DD experiments. Our work expands the existing DM direct detection literature in two important ways: a) we consider the effect of astrophysical uncertainties for DM models that encompass a diverse set of operators in non-relativistic effective field theory (NREFT) with different recoil energy and velocity dependences; b) we adopt a novel approach to analyze the DM phenomenology at DD experiments through just two experimental observables: the total number of events and the number of events per bin, or the shape, of the recoil spectrum. These ingredients are then self-consistently combined in our statistical framework that uses the euclideanized signal (ES) method developed by refs. [40,41] to forecast the ability of next-generation experiments to resolve DM model parameters using SHM and Gaia velocity distributions as two representative cases. 1 We note that our method is quite different from both traditional forecasting approaches for DM direct detection that relied on expensive Monte Carlo (MC) simulations for generating mock data sets [42][43][44][45][46], or analyses which studied the effect of uncertainties in SHM on constraining the DM particle physics properties [47,48].
We show that the Gaia velocity distribution could significantly change not just the overall rate, but also the shape of the DM recoil spectrum at DD experiments. Moreover, we find that both the DM mass and the recoil energy dependence of the model could enhance or suppress the effect of the Gaia velocity distribution vis-á-vis the SHM while inferring the DM model parameters. This is another example of the curious interplay between different DM properties illustrated in Fig. 1 that, to the best of our knowledge, has not been fully explored. More concretely, we study the effects of DM velocity distribution on: i) a simultaneous reconstruction of DM mass and coupling, as well as ii) mediator and DM masses, for a given model, and iii) prospect of pairwise model discrimination. We also revisit the idea of combining DD data from experiments with different targets to reduce the effect of astrophysical uncertainties and, in some cases, dramatically improve sensitivity to certain parts of DM mass-coupling parameter space.
The paper is organized as follows. We review the basic ingredients to compute the recoil spectra at DD in section II: we present the main equation for differential recoil rate in section II A, summarize and review the velocity distribution taking into account of the Gaia Sausage in section II B, and review the NREFT formalism along with several benchmark DM models in section II C. In section III, we demonstrate how Gaia velocity distribution could modify both the overall rate and the shape of the recoil spectra in the NREFT framework. In Section IV, we introduce the recently developed ES statistical framework, which allows us to make forecasts without running MC simulations. In section V, we present our results on the effects of Gaia distribution on reconstructing DM model parameters and distinguishing between different models at future DD experiments as well as setting constraints using current data. We conclude in section VI.

II. PHENOMENOLOGY OF DARK MATTER-NUCLEAR INTERACTION
In this section, we briefly review the key ingredients to compute the rate of DM scattering off nucleus in a DD experiment. In particular, we will discuss possible new DM velocity distributions inferred from the Gaia data. We will also review both the model-independent framework, the non-relativistic effective theory and some specific benchmark models to study different types of DM-nucleus interactions.

A. DD basics
The rate of DM scattering off nucleus in a DD experiment is given by R = n χ σv χ , where n χ and v χ are the local DM number density and speed relative to the Earth respectively, σ is the scattering cross section and . . . indicates an average over the local DM velocity distribution. The differential recoil rate for a target T per unit recoiling energy can be written as, where in the second line, we use and group together different factors by the main type of physics they rely on. In the master equation for the differential rate, ξ T is the mass fraction for each type of target nucleus T (the detector could be composed of different nuclides), m T is the target nucleus mass and N T = 1/m T is number of scattering centers per unit mass. ρ χ is the local DM density in the solar system, which we take to be ρ χ = 0.4 GeV/cm 3 [49][50][51]. 2f ( v) is the local DM velocity distribution in the Earth frame. In the next section, we discuss in some detail different possible velocity distributions, both the SHM, which has been most commonly used, and the new empirical distributions constructed from Gaia data [31,33,35,52]. The velocity integration range is bounded from above by the escape velocity of DM particles v esc . The minimal velocity for DM to scatter with recoiling energy E R , in the case of elastic scattering, is where µ T is the DM-nucleus reduced mass. We only discuss elastic scattering throughout the paper. Lastly, in the part that depends on particle and nuclear physics, J and J χ are nuclear and DM spins respectively. |M| 2 is the scattering matrix element squared averaged over 2J χ + 1 and 2J + 1 initial DM and nuclear spins, and summed over the final spins. We postpone a more detailed discussion of the DM-nuclear scattering theory to Sec. II C. Note that the recoiling spectrum depends on the detector material in multiple ways. For example, in the velocity integration, v min depends on the target nucleus mass and the matrix element depends on the nuclear form factor, which is determined by the type of the target.

B. DM velocity distribution
The local velocity distribution in the GC frame for the SHM, i.e. for an isotropic, isothermal DM halo in equilibrium, is well-modeled by a truncated Maxwell-Boltzmann distribution [8], where we take the velocity dispersion and escape velocity of the DM halo to be σ v ≈ 160 km/s and v esc ≈ 540 km/s respectively [53], and the normalization constant N esc is given by, Since we are interested in estimating the DM recoil rate for a detector on Earth, we need to boost the DM velocity distribution by the Earth's velocity in the galactic frame, v obs (t), At the leading order, the DM recoil spectrum can be decomposed into a constant and a periodic component. The constant component arises from the velocity of DM halo in the heliocentric frame, whereas the Earth's orbital motion around the Sun manifests as a sinusoidally varying component, commonly referred to as annual modulation [54,55]. We ignore annual modulation effects in our analysis, assuming a time averaged value for the Earth's velocity v obs ≈ 230 km/s instead. A crucial assumption in the formulation of the SHM is the condition of local equilibrium. However, if the MW has undergone one or more recent mergers, the equilibrium condition is then invalid, and we would need a method to empirically determine the DM velocity distribution. Assuming the CDM paradigm, hierarchical structure formation implies that the MW halo should primarily consist of virialized tidal debris from old subhalo mergers with other spatial and kinematic substructure sourced by more recent ones. Ref. [56] used the Eris N -body simulation to show the correlation between velocities of old, metal-poor stars and the virialized DM component of a MW-like halo. Along the same lines, ref. [15] argued that the velocity distribution of stars in a class of substructure called debris flow 3 is a good kinematic tracer of its accreted DM counterpart based on the Via Lactea simulation.
In unarguably the golden-age of data-driven astrophysics, we can now obtain 7D information 4 for main-sequence stars in the MW halo by cross-matching the Gaia data releases [25,26] with the Sloan Digital Sky Survey (SDSS). Using this cross-matched catalog, refs. [27,28] found signatures of a debris flow in the solar neighborhood (within ∼ 4 kpc of the Sun)-the so-called Gaia-Sausage -that consists of metal-rich halo stars with high radial anisotropy.
Further investigation of other phase space substructures [27,29,33] indicates that, in fact, the MW might have experienced at least two different accretion events, namely those of the Sausage [27] and Sequoia [58] dwarf galaxies. The Gaia-Enceladus structure [28] hints at a possible third event, although it appears to partially consist of debris from the other two mergers [52].
For our analysis, we only focus on the effects of the Gaia-Sausage, since it is the dominant merger in the solar neighborhood contributing ∼ 70% of all accreted stars [27]. Ref. [31] used the SDSS-Gaia DR2 data set for a subsample of ∼ 190, 000 stars to make the first empirical determination of the local DM velocity distribution. They performed a Gaussian mixture model (GMM) analysis on the joint distribution of stellar velocities and metallicities to classify stars in three populations with distinct kinematic properties: metal-rich young disk stars formed in-situ; accreted stars which include metal-poor stars in the smooth isotropic halo, and intermediate metallicity stars with a high radial anisotropy that constitute the Gaia-Sausage substructure. The velocity distribution of each component is then determined by sampling the joint distribution assuming broad, uniform priors on the GMM parameters.
The total DM velocity can then be written as a linear combination of the substructure and halo velocities, weighted by the the ratio of DM in each component. In terms of velocity distributions, where η sub parametrizes the DM fraction in substructure. The fraction of stars in substructure relative to total accreted stars is typically a poor tracer for η sub as simulations show that the halo DM component may contain significant contributions from accretion of diffuse DM and DM in non-luminous subhalos. In fact, as demonstrated by ref. [32,59], late-time accretion from the latter component may affect both the shape and the velocity distribution of all accreted DM. We defer a detailed study of these effects to future work, and adopt a distribution for η sub (shown in the left panel of Fig. 2) estimated in ref. [32] using an approximate relation between the mass-to-light ratio and metallicity for the MW [60,61]. We plot the Gaia (shorthand for empirical) DM velocity distribution in the heliocentric frame for each component 5 in the right panel of Fig. 2 assuming a median value of the DM substructure fraction, η sub = 0.42. We emphasize that there is no smooth interpolation in η sub between the SHM and the Gaia velocity distribution, i.e. η sub = 0 does not yield the SHM. Finally, we comment on ref. [63], a very interesting study of the high-velocity DM particles that appeared recently. The authors found, using an idealized N -body simulation, that the DM halo of the Large Magellanic Cloud (LMC) could both accelerate the local MW DM particles to higher velocities and contribute its own high velocity particles. This would result in the local DM velocity distribution to be skewed towards high velocities as compared to the SHM, thereby strengthening cross-section limits at low DM masses. While we find their premise intriguing, studying the effect of LMC on DM in the Sausage debris flow requires a more thorough investigation with cosmological zoom-in simulations, which is beyond the scope of our analysis.

C. DM-nucleus scattering theory
In this section, we will briefly review a model-independent framework and some concrete models to study different types of DM scattering off nucleus in DD experiments. There is a huge literature on possible DM scattering in DD and we do not intend to provide an exhaustive review here. We will only refer the reader to the original papers and papers we actually use for our analysis.
1. Model-independent framework: non-relativistic effective field theory The typical DM velocity in our galaxy is v ∼ < 10 −3 c. The incident DM kinetic energy is around O(10 keV) for a DM particle with mass ∼ 10 GeV. Given that DM scattering in DD experiments is non-relativistic (NR), a simple model-independent way to parametrize different types of DM-nucleon interaction is the Galilean-invariant NR effective field theory (NREFT), first proposed and developed in refs. [64] and [65]. The core result is that in the NR limit, DMnucleon interactions could be encoded in 16 NREFT operators, 15 of which are linearly independent. These operators are expressed in terms of four three-vectors: DM spin S χ , nuclear spin S N , the momentum transfer q = p − p with p ( p ) the incoming (outgoing) DM three-momentum and the transverse velocity We will use only 12 linearly independent operators, O i 's, listed in Table III of Appendix A. These 12 operators are usually sufficient to describe the NR limit of many relativistic operators that appear in simple models with spin-0 or spin-1 mediators. In the NR limit, a relativistic operator in the field theory can be mapped onto a linear combination of the NR operators. Thus a relativistic Lagrangian for a particular DM scattering model, which could contain several relativistic operators for DM-nucleon interaction, could be written in terms of the NREFT operators as, where N = n, p labels the type of nucleon DM interacts with, which could be either neutron or proton and i labels NREFT operators. The coefficients, c i 's, depend on the coupling coefficients in the relativistic theory and the Wilson coefficients obtained by mapping the relativistic operators to the NR ones. Compared to the standard literature (e.g., refs [65][66][67]), we also take the mediator propagator out of c i 's and explicitly write it out (m med;i is the mediator mass for the ith interaction).
In the paper, we only consider the leading order spin-independent (SI) scattering in which DM scatters off the entire nucleus coherently. To go from DM scattering with individual nucleons to scattering with nucleus, one needs to take into account of the nuclear response which is encoded in the form factor Further development of effective field theory from quarks to nucleons could be found in refs [68,69]. The spin-averaged amplitude squared is then Note that our normalization of the form factors F's differs from that in the literature (e.g., refs [65,66]) by a factor of (4m χ m N ) 2 with m N the nucleon mass. In NREFT, the differential scattering rate in Eq. (1), in terms of the c i 's and the form factors, is The form factors used in our analysis are provided in Appendix A.

Some benchmark models
A DM model could contain several different DM-nucleon interactions and a relativistic operator between DM and nucleons could map onto multiple NR operators. Interpretation of the data for a specific DM model may not be a simple combination of that for each individual NR operator it could map onto. For instance, the constraint on the coupling of the relativistic operator could be a complicated sum of constraints on each NR operator it corresponds to with different weights due to interferences between different NR operators. Thus in analyzing how the new velocity distribution affects interpretation of DD data, we also consider a few simple representative models. Again we do not intend to be exhaustive. We only select and compare a few models with different E R and velocity dependences, which could be affected by SHM and Gaia velocity distributions in different ways.
We assume for simplicity that the DM is fermionic. We analyze six SI DM models: DM interacting through heavy scalar mediator (leading to the simplest contact interaction), millicharged DM (mC), DM with magnetic dipole moment (MD) with either heavy or light mediators, DM with electric dipole moment (ED) with light mediator and anapole DM with a heavy mediator. All the models are summarized in table I.

Model
Relativistic Operator NREFT Operator ER and DM velocity moment Contact interaction (heavy scalar mediator) gcχχN N gcO Magnetic dipole (MD) w/ heavy mediator We also show the leading order E R and velocity moment dependences of the corresponding spectra. The coupling constant in front of the operator defines the model parameter to be constrained, e.g. coupling strength g c for contact interaction, charge fraction χ for millicharged DM, and the dipole moment µ χ for magnetic dipole DM.
Heavy scalar mediator: the model is simply where φ is the heavy scalar that could be integrated out, giving rise to a four fermion operatorχχN N . One classic example is the Higgs-portal model in which the heavy scalar mediator is the SM Higgs particle [70]. Millicharged DM: DM carries a small electric charge χ e and couples to the SM photon, A µ , χ particle can either carry a fractional electric charge directly [71][72][73] or charged under other U (1) which kinetically mixes with the SM hypercharge [74]. A µ couples to the SM current, which is given by where g p and g n are the g-factors of protons and neutrons respectively. Magnetic dipole DM with a light mediator: DM does not carry electric charge, but it could still couple to the SM photon through a loop of other charged species. This could generate an anomalous magnetic dipole moment of the DM particle [75]: Magnetic dipole DM with a heavy mediator: magnetic dipole DM can also interact with the SM electromagnetic current through a heavy mediator. One way is through the kinetic mixing of a broken dark where the dark photon is heavy with a mass larger than the momentum transfer in DD, m 2 D q 2 . Integrating out the heavy dark photon leads to a coupling between DM and the SM photon as listed in Table I. Electric dipole DM with a light mediator: if a dark sector has CP violating interaction, it is also possible that the DM carries an electric dipole moment which aligns with the DM spin. Anapole DM with a heavy mediator: if DM is a majorana fermion, monopole and dipole interactions with electromagnetism is forbidden by CPT symmetry. The only allowed EM coupling is the anapole coupling. For example, the majorana DM could have the anapole moment under a dark broken U (1), which kinetically mixes with the SM photon [76]. Integrating out the heavy dark photon, we have Early references on anapole DM could be found in [77,78].

III. EFFECT OF THE DM VELOCITY DISTRIBUTION ON DD RECOIL SPECTRUM
Before presenting the statistical method for DD forecasting and the results, we want to discuss how different DM velocity distributions could affect the differential recoil spectrum, the key quantity in DD experiments. This will help us obtain some qualitative ideas and physical intuition of the effects of DM velocity distribution on interpreting current and future DD data.
The form factor F i,j could depend on the incident DM velocity squared v 2 (as well as momentum transfer squared q 2 ) for DM-nucleon interaction beyond the simplest contact one, which could contribute to the power of velocity in the velocity integration. Schematically the form factor could be written as a power series of v 2 , where F (n) (q 2 ) is the nth q-dependent coefficient in the Taylor series of v 2 . Plugging Eq. (19) into Eq. (12), we have where While higher order velocity moments could be present, we only need to consider n = 0, 1 for most DM models in the literature. Note that these velocity integrals are the partial moments of f (v) given that the integration starts from v min as opposed to from 0 and the upper end of the integration is bounded by v esc instead of ∞. In the computation presented below, we will take v esc to be ∞ for the Gaia velocity distribution to be consistent with ref. [31]. We present g(v min ) and h(v min ) as functions of v min in Fig. 3. No matter what the velocity distribution is, g(v min ) and h(v min ) are always monotonic decreasing functions of v min . The differences between the ones derived from different velocity distributions originate from the fact that the Gaia velocity distribution peaks at a lower v compared to the SHM distribution and thus has a smaller tail at high velocity, as one could see in the right panel of Fig. 2. For g(v min ), the weight in the velocity integration is 1/v, which also prefers small velocity. Thus for small v min , g(v min ) derived using f Gaia is larger than that from f SHM . Yet when v min increases to 250 km/s, the smaller high velocity tail of f Gaia becomes the determinant factor that makes g(v min ) Gaia smaller than g(v min ) SHM . For h(v min ), the weight in the velocity integration is v, which prefers large velocity. In this case, even for small v min , the weight overcomes f Gaia 's preference for small velocity. At larger v min , both the weight and the smaller high velocity tail of f Gaia suppresses has been multiplied (divided) by the speed of light to be made dimensionless, for a better illustration of their relative magnitudes.
The recoil spectrum depends on not only the velocity moments but also the coefficients F i,j (q 2 ) is, in general, a polynomial of q 2 multiplying an exponent suppression factor ∼ e −q 2 s 2 , in which the effective nuclear radius is s ≈ 0.9 fm. Since the exponential factor is common for all nuclei, we will focus on the polynomial part of F (n) i,j (q 2 ) that depends on the type of DM-nucleon interaction. The momentum transfer is related to the recoiling energy E R and v min as, We could express q 2 in terms of v min and the differential recoil rate as a function of v min only. To illustrate how the recoil spectra vary with velocity distributions, we select three NR operators as examples and compute their recoil spectra assuming a xenon target. For simplicity, we consider only one NR operator at a time with a single type of mediator. Let's consider two limits to illustrate the effects of different DM velocity distributions: heavy mediator with m 2 med q 2 so that the propagator is approximately a constant, and light mediator with m 2 med q 2 so that the propagator is approximately 1/q 2 . It is not difficult to generalize the discussion to cases in between the two limits with m 2 med ∼ q 2 . In the two limits, we have the dependence of the recoil spectrum on v min for each operator schematically as To derive the spectrum's dependence on v min for O 8 , we use v ⊥ = v + q 2µ T and v ⊥ · q = 0. We could understand some general features of the spectra independent of the velocity distributions. Given the exponential factor in the nuclear response function, all spectra fall off at large recoiling energies. In the heavy In this figure and the next one, the differential recoiling rate is in arbitrary unit.
mediator limit, the spectrum of O 1 peaks at low recoil energy (or equivalently, small v min ) since g(v min ) is a monotonic decreasing function of v min . The differential rate of O 11 contains an additional factor, v 2 min , which prefers larger E R . Thus the spectrum peaks at a higher E R , away from the detection threshold, due to a balance between v 2 min and g(v min ). For O 8 , the differential rate is a sum of two terms with opposite behavior: v 2 min g(v min ) peaks at large recoiling energy while h(v min ) peaks at threshold. Numerically it turns out that h(v min ) dominates over v 2 min g(v min ) so that the spectrum still peaks at low recoil energy. The situation is much simpler in the light mediator limit. The spectra for all operators peak at threshold due to the mediator propagator 1/q 2 .
We present the spectra for both light and heavy DM with masses at 12 and 50 GeV respectively in Fig. 4 (in the heavy mediator limit) and Fig. 5 (in the light mediator limit). In the plots, the cyan curves are based on SHM while the orange curves are based on the Gaia distribution with η sub = 0.42. Based on the discussion above, we could understand further some details of the spectra for each velocity distribution.
In the heavy mediator case, • For both O 1 and O 11 , spectra based on Gaia peaks at a higher value at low recoiling energy and falls off faster compared to the one from SHM given the shapes of g(v min ) shown in Fig. 3a.
• For O 1 , the spectrum based on Gaia distribution is steeper for light DM than that for heavy DM. More specifically, the recoiling energy at which the differential rate from SHM becomes larger than that from Gaia for heavy DM (the cross-over of the two curves) is larger than that for light DM. This is because at v min ∼ 250 km/s, the relative sizes of g(v min )'s for Gaia distribution and SHM switches (Fig. 3a). For a given v min , the more heavy DM is, the larger E R it corresponds to, as one could see from Eq. (23). Similar argument could be used to explain the differences for spectral of light and heavy DM scattering through O 11 .
• For O 8 , the spectra based on Gaia are always below those based on SHM since the scattering is mostly determined by h(v min ).
Similar features are present in the light mediator case (though they are less evident in Fig. 5): • For O 1 and O 11 , the Gaia spectra are steeper than those of SHM. At low E R , the differential rate based on the Gaia distribution is greater than that based on SHM while the opposite is true at higher E R . Analogous to the heavy mediator case, the cross-over happens at a larger E R for heavy DM compared to light DM.
• For O 8 , Gaia's differential rate is always below the SHM's mainly due to the dominance of h(v min ).
In general, the overall rate of light DM with mass 10 GeV could be suppressed with a Gaia distribution compared to that of SHM. This leads to weaker constraints and poorer determinations of the light DM parameters using Gaia velocity distribution. On the other hand, for heavy DM, the relative sizes of scattering rate with either Gaia or SHM distributions depends on the type of the interaction. The scattering rate with Gaia could be enhanced when the associated velocity moment is g(v min ) and the scattering rate is proportional to a non-negative power of v min (or equivalently, q). This could lead to stronger constraints and better determinations of the DM parameters, e.g., when the model maps onto O 1 and O 11 with heavy mediators. For interactions associated with h(v min ), however, the recoil rates assuming Gaia distribution are suppressed compared to that of SHM, resulting in a weaker constraint on the coupling. Lastly, the differences of the spectral shapes in the light mediator case, which depend on negative powers of v min (q), with either Gaia or SHM distributions are small. These qualitative discussions based on the recoil spectra will indeed be confirmed with numerical computations in Sec. V.

IV. STATISTICAL FRAMEWORK
In this section, we introduce the statistical framework to study the effect of uncertainties in the DM velocity distribution while reconstructing particle physics parameters with DD data. As discussed previously in the introduction, for a given DD experiment, there is a three-fold degeneracy between the different classes of DM parameters. At the same time, from a statistical inference perspective, we can only access (a combination of) these parameters through experimental observables. For DM-nuclear interactions, these observables are simply the overall rate and the number of events per bin or the shape of the recoil spectrum (which may also include background events). The qualitative relationships between these observables and DM parameters of interest are summarized in table II.  However, any statistically significant number of DM events have not been detected at Pico [79], Lux [80], Super-CDMS [81], PandaX-II [82], and Xenon-1T [83]. Thus, for a given data set, recent experiments such as Xenon-1T usually report the number of events for a best-fit WIMP candidate at some reference mass and cross-section assuming contact interaction with a heavy mediator, alongside information regarding the backgrounds and experimental efficiency [53,83]. 7 This is the main reason why most recent theoretical analyses that have investigated the effect of an empirical DM velocity distribution on constraining its particle physics properties [31,[36][37][38]48] use marginalized upper limits -a method that relies on comparing the total number of predicted and observed events and doesn't include any information about the spectral shape -to indicate their results.
We circumvent the restrictions that come from using existing experimental data by forecasting the sensitivity of nextgeneration experiments to DM model parameters in presence of astrophysical uncertainties. Traditional forecasting approaches for DM direct detection leverage both the overall rate and shape spectrum for recoil spectra by generating mock data sets with Monte-Carlo (MC) simulations. The degeneracies in model parameters are then extracted by either calculating the profile likelihood [42,47], or sampling over their posterior distribution [43][44][45].
The major drawback of using MC simulations is that they are computationally expensive and rely on the choice of several pre-defined benchmark points. In fact, a natural way of forecasting the degeneracies between DM model parameters is to simply compare -in a statistically consistent way -the spectra from various points in the parameter space. Along these lines, refs. [40,41] introduced a novel method based on information theoretic techniques that facilitates fast, benchmark-free model comparison. A python implementation of their results is available in the open source swordfish code, 8 and a proof-of-concept application to future DD searches has been demonstrated in ref. [86]. In the rest of section, we summarize the main idea behind the euclideanized signal (ES) method followed by a short discussion on how to interpret its results with a concrete example.
Consider a d-dimensional parameter space, θ ∈ Ω M ∈ R d for a given DM model M, and the associated Asimov data set D A (θ) for each point [87]. 9 Given two model parameter points θ 1,2 ∈ Ω M , we can construct a likelihood 7 We realize that this is a rather simplistic way of describing the highly non-trivial procedures used by experimental collaborations in analyzing their data. Our main point is that, short of performing the full event-by-event likelihood analysis described in refs. [84,85], we as theorists currently only have access to the total number of events for the DM signal spectrum but not its shape. 8 Code for [41]: https://github.com/cweniger/swordfish/blob/master/Examples/swordfish_DD.ipynb 9 We emphasize that there is no mock data generated at any stage of our analysis. The use of Asimov data set in eq. (27) implies that the sampling distribution gives the median significance for two hypothetical data sets which have the parameter points θ 1,2 as their maximum likelihood estimates. ratio test statistic (TS) [88] with well-defined asymptotic behavior [89], which we use to reject the null hypothesis that θ 1 and θ 2 are indistinguishable at the (1 − α)% confidence level (CL).
I is the profiled Fisher information matrix, which we explicitly derive in Appendix B. The summation runs from 1 to n b , the number of data bins. Commonly referred to as Fisher-Matrix forecasting, this approach is extensively used for making sensitivity projections of a model experiment in cosmology and astrophysics. Unfortunately, calculating the TS in eq. (26) for model comparison of N points in the parameter space can be computationally expensive when N is large, since it requires O(N 2 ) matrix inversions. A further complication arises due to the presence of derivative terms in the calculation ofĨ that could introduce additional numerical noise. Noting that the profiled Fisher information I transforms as a metric on the parameter space θ, ref. [41] was able to reduce the number of matrix inversions to O(N ) by mapping the parameter space into a higher-dimensional signal space and expressing the TS as a euclidean distance between two signals. More concretely, the θ → x( θ) ∈ R n b , transforms the parameter space to the n b -dimensional signal space with unit Fisher information matrix. Details of the embedding can be found in the Appendix B. After this transformation we can express the TS in terms of the euclidean signal x i , Eq. (27) is the main ingredient of our benchmark-free forcasting approach, and from where the ES method obtains its name. In the language of this method, as long as the parameter space is sufficiently sampled, signal discrimination is only possible at the (1 − α)% CL if the signals from two parameter points are at least a distance r α (M) apart in the projected signal space. The distance, in turn, is related to the sampling distribution of the TS, where F χ 2 d is the cumulative distribution function (CDF) of the χ 2 distribution. The above procedure can be understood very loosely as comparing two distributions, albeit incorporating the fact that they arise from the same likelihood function. We illustrate the efficacy of this method in Fig. 6, the left panel of which shows the constraints in DM mass-coupling space for a model with contact interaction mediated by a heavy scalar particle. The closed ellipses represent the usual 68% CL contours in parameter space for arbitrary benchmark points at a next-generation Darwin-like liquid Xenon (LXe) experiment [6] assuming SHM. These are obtained by constructing hyperspheres of radius, r α , in the Euclidean signal space, and back-projecting them to the parameter space using a lookup table for the embedding map. For a χ 2 d distribution with d = 2, eq. (28) implies that 68% CL corresponds to a threshold value of r 0.32 = 1.52.

V. RESULTS
In this section, we put together the formulas and methodology presented in previous sections and use them to study how DM velocity distribution inferred from Gaia sausage could affect reconstruction of various DM particle physics parameters at next-generation DD experiments. For concreteness, we only consider a Darwin-like liquid Xenon (LXe) experiment and a complementary DarkSide-20k-like Argon experiment with both high [7] and low mass [90] search programs (see Appendix C for more details). Unlike Sec. III, we present our results for the benchmark DM models listed in Table I instead of individual operator in NREFT. While examining recoil spectrum of each operator is insightful, concrete models, especially those with well-motivated UV completions and/or distinct phenomenologies, enable an easy comparison of our results with those in the literature. In addition, since there could be non-trivial mapping between a model and NREFT, it may not be straightforward to find the sensitivity of DD to a full-fledged model by combining the sensitivity to individual NR operators. Yet we will still find the qualitative understanding developed in Sec. III a useful starting point for the results discussed here.

A. DM Mass -coupling
We forecast the sensitivity of a next-generation LXe, Darwin-like, experiment to simultaneously resolve the DM mass and coupling for SHM and Gaia velocity distributions. Our first step is to investigate the effect of uncertainty in  the DM substructure fraction η sub , estimated by ref. [32] to be η sub = 0.42 +0.26 −0.22 . We show, in the right panel of Fig. 6, the 68% CL contours and 90% CL upper limits corresponding to a Gaia velocity distribution with the median and ±1σ DM substructure fractions. Comparing them with constraints for the SHM leads us to conclude: the primary effect in reconstructing DM model parameters arises due to the qualitative differences between the SHM and Gaia velocity distributions, while the variation of the substructure fraction is only a subdominant effect beyond it. Thus, for the rest of our analysis, we fix the DM substructure fraction to its median value, η sub = 0.42. 10 In Fig. 7, we show the 68% CL contours obtained using the ES method for four DM models with characteristic E R and velocity moment dependences (outlined in table I), namely contact interaction, millicharged DM, and DM with electric and magnetic dipole moments, for SHM and Gaia velocity distributions. For reference, we also plot the 90% CL upper limits following the latest Xenon-1T results [83] and projected upper limits for future LXe experiments using the equivalent counts method [40,41]. We constrain couplings instead of cross-sections (cf. [45,86]) as we don't integrate over the entire E R range to obtain the respective cross-sections for the DM models we consider here. Although, as illustrated in the left panel of Fig. 6, the ES method allows us to plot the degeneracy contour for any point in the parameter space, we show our results at two benchmark points corresponding to light (m χ = 12 GeV) and heavy (m χ = 50 GeV) DM for easier interpretability of our results. To ensure we are making an apples-to-apples comparison when studying the changes in constraints across models, we choose couplings such that the number of events is the same for each benchmark point with SHM.
Before discussing the effect of DM velocity distribution, we explain the general behavior of constraints in DM coupling-mass space given in Fig. 7 in terms of the recoil spectra shape and the total event rate. For a given DM model, the recoil spectra for low mass DM peaks closer to threshold than for heavy DM. Moreover, the shape of the recoil spectra is degenerate only for a narrow range of masses, whereas a change in the total rate can be compensated by a wide range of couplings. As a result, in the light DM regime, we observe a large degeneracy in the coupling but a reasonable DM mass resolution. Conversely, for heavier DM, a DD experiment is more sensitive to the couplings but suffers from poor mass resolution, because the recoil spectra shape is degenerate for a wide range of DM mass. Thus, given a similar number of events, we can heuristically treat the light and heavy DM regimes as shape-and total rate-limited respectively. While these interpretations hold generally, certain qualitative details like the shapes and sizes of the contours can vary significantly between the SHM and Gaia velocity distributions. For instance, in the light DM regime, we observe that the contour size increases for all models with Gaia velocities, implying that LXe experiments have reduced sensitivity to Gaia distribution as compared to the SHM. This observation is consistent with the weakening of upper limits on SI cross-section for light DM first reported by ref. [31] with the latest Xenon-1T data. Our results also indicate an interesting effect of Gaia velocity distribution that has not been previously discussed in the literature: depending on the model, there is a marginal improvement in the sensitivity of LXe experiments to heavy DM for Gaia when contrasted with the SHM. We also note that for models with light mediators, the experimental sensitivity becomes poorer across the entire mass range irrespective of the DM velocity distribution.
These results can be understood, at least to leading order, in terms of the E R and velocity moment dependences of each model. We start with the observation that contours of the Gaia velocity distribution are less constraining than those of SHM at low DM masses for all models we consider here. These models have a leading order DM velocity moment that scales as g(v min ) and/or h(v min ) suppressed or enhanced by additional powers of E R . Since light DM corresponds to a high v min for a heavy target like xenon, only the tail of the g(v min ) and/or h(v min ) distribution (Fig. 3) contributes to the recoil rate, where the SHM curves always dominate over the Gaia ones.
For heavy DM, on the other hand, varying DM mass could lead to sharper changes in the recoil spectra shapes with the Gaia distribution as compared to the SHM. Thus, there is an improvement in the sensitivity to DM mass as evidenced by the shrinking of 68% CL forecast contours in the mass direction in Fig. 7. Moreover, as shown in the top row, this effect is most apparent for models with a heavy mediator, or equivalently for non-negative powers of E R . We also note a subtle difference between models with contact interaction and magnetic dipole interaction. For contact interaction, the upper limits for DM with mass above 50 GeV are slightly tightened with the Gaia distribution. This is due to an enhanced recoil rate contributed by Gaia's larger g(v min ) at v min ∼150-200 km/s as compared to SHM. Whereas, in the case of magnetic dipole DM, the Gaia upper limits for DM with mass above 50 GeV are slightly weakened, since the scattering rate of magnetic dipole DM scales as E 2 R g(v min ) + E R h(v min ). The second term proportional to h(v min ) leads to a small reduction in the overall recoil rate with the Gaia distribution, as compared to the SHM one. Meanwhile, the positive powers of E R for magnetic dipole interaction result in an enhanced sensitivity of forecasts using the Gaia distribution. This is an interesting example where despite the reduction in the total number of events, the sensitivity actually improves with the Gaia distribution! In case of light mediators, the inverse powers of E R make the recoil spectra peak sharply as E R → 0. Yet for a finite threshold E R ≈ 5 keV, the Darwin-like experiment is only sensitive to the tail resulting in highly degenerate recoil spectra for different velocity distributions. This leads to poorer experimental sensitivity for both SHM and Gaia velocity distributions across all DM masses for models with light mediators, compared to the contact interaction model with a heavy mediator. In this case, lowering the detection threshold could improve the sensitivity to light mediators, making another physical case for the low threshold frontier.
Notice that we have listed six models in Sec. II C 2 and we only show results for four models in Fig. 7. Results for the other two models, anapole DM with a heavy mediator and magnetic dipole DM with a light mediator are presented in Appendix D. They have very similar features to the results in this section and it is straightforward to understand those following the discussion above. DD experiments are best suited to constrain the mediator mass m med when it is at the same order of the momentum exchange, i.e. m med ∼ q ∼ O(10 MeV) [47,91,92]. We use the ES method to simultaneously constrain the mediator and DM masses at a Darwin-like experiment for a fixed coupling. We also study the effect of the DM velocity distribution on the forecast. For simplicity, we only consider DM with contact interaction mediated by a light scalar particle and fix the DM coupling to be the same for all benchmark points.
In the left panel of Fig. 8, we show the 68% CL forecast contours for SHM in the m −1 med − m χ plane. Broadly, the structure of these contours resembles the fishnet plot in the DM mass-coupling plane shown in the left panel of Fig. 6. On closer inspection, however, we can roughly delineate three regimes of sensitivity in this parameter space. The upper part of the plots is the light mediator (m med q) regime where the propagator squared simply scales as ∼1/q 4 . In this case, the m med dependence drops out, and DD experiment is insensitive to the mediator mass. Next, we consider the heavy DM regime in the lower right part of the plots, where we have chosen the mediator mass such that m med q. We find that the degeneracy in the recoil spectra in this limit is due to the DM mass, and any change in the mediator mass effectively acts as a rescaling of the overall coupling. Lastly, for intermediate DM mass and m med q, our Darwin-like experiment can precisely reconstruct both mediator and DM masses primarily due to the high number of signal events (a factor of a few greater than the other benchmark points) in this regime. where the light-colored bands indicate the 1σ Poisson uncertainties. The third column shows the 68% CL forecast contours for SHM (cyan) and Gaia (orange) in the coupling-coupling space for two comparison models: (first row) millicharge with a light mediator and contact interaction with a heavy mediator, (second row) millicharge and magnetic dipole both with a light mediator (SM photon). We also include 90% CL upper limits from the latest Xenon-1T results (yellow) and projected upper limits for a Darwin-like experiment (indigo) assuming SHM (solid) and Gaia (dashed) velocity distributions.
We also illustrate the differences between SHM and Gaia velocity distributions in the right panel of Fig. 8  In presence of a positive signal at a future DD experiment, one of the most important goals is to determine the type of DM-nuclear interaction and discriminate between different model candidates. To demonstrate the model selection, we postulate a scenario in which there are two candidate models of interest. We parameterize our model as the sum of a pair of interactions, While holding m χ and m med fixed, we sample different values of (c a , c b ) and test how well a given recoil spectrum shape can determine the model parameters at 68% confidence level. We test the two pairs of models: 1) a, b = contact interaction, millicharge, and 2) a, b = magnetic DM with light mediator and millicharge respectively. We present the results in Fig. 9 for light DM with mass at 12 GeV and Fig. 10 for heavy DM with mass at 50 GeV. These results are all based on a Darwin-like LXe experiment. In Fig. 9, we present 68% CL forecast contours in the coupling-coupling space for a 12 GeV DM particle. In the top row, we have a millicharged DM giving rise to an experimental signal and we want to test whether we would confuse it with the simplest contact interaction as both interactions lead to spectra peaking at experimental threshold. This scenario is equivalent to setting c a = 0, c b = 0 in eq. (29), where c a = χ e for the millicharge model and c b = g c for the contact interaction. In the bottom row, we have a DM particle with both magnetic dipole moment and millicharge interacting with the nucleus through the SM photon contributing with comparable rates, where c a = eµ χ /2 and c b = χ e. We want to test how well we could constrain the two relevant electromagnetic moments.
From the third column, one could see that while it is possible to reconstruct the model parameters (with large uncertainties) assuming SHM, the discrimination power is entirely lost with the Gaia distribution. The millicharged light DM could be misidentified as a light DM with simple contact interaction at 68% C.L, as shown in the last plot in the top row. DM with both millicharge and magnetic dipole moment could not be distinguished from DM with only one of them, as shown in the last plot in the bottom row. This result could be understood from the recoil spectra shown in the first two columns, in which we fix the couplings of different DM models to give the same event numbers with SHM. Comparing the spectra based on the Gaia velocity distribution with those from SHM, one find that with Gaia distribution: i) the total number of events is significantly lower, which, in turn, increases the Poisson uncertainty; ii) the spectral shapes of different models, especially the tails of the distributions, are more degenerate. This is consistent with what we find in Sec. III using NREFT.
Results of model discrimination for heavy DM with mass at 50 GeV are presented in Fig. 10. In this case, we find that reconstructing the couplings is more accurate as compared to the light DM case, independent of the model combinations. We could almost distinguish different models or different combinations of model parameters equally well for both SHM and the Gaia distribution. The main reason is that the total number of events and the spectral shape do not change much when the velocity distribution varies. The idea of combining different targets for a more accurate identification of momentum dependence of DM interactions, or for precise reconstruction of DM mass and model parameters is well-documented in the literature [42,45,[93][94][95]. We revisit this idea by using the ES method to forecast the sensitivity of two complementary next-generation experiments with different targets. In particular, we focus on a Darwin-like and a DarkSide-20k-like experiment with xenon and argon targets respectively. Besides forecasting for high mass DM as officially proposed by the DarkSide-20k collaboration [7], we also include results for a dedicated low DM mass search with a low threshold configuration similar to ref. [90].
For simplicity, we consider the contact interaction as an example. Fig. 11 shows the 68% CL contours in the DM mass-coupling space for each of these experiments at two benchmark points corresponding to light (m χ = 12 GeV) and heavy (m χ = 50 GeV) DM. We find that combining forecasts from both xenon and argon targets could dramatically improve the sensitivity for low mass DM and, to a lesser extent, high mass DM. This result is independent of the model for DM velocity distribution. The extremely high resolution for low mass DM is due to the low threshold version of a DarkSide-20k-like experiment as listed in table IV of Appendix C. An advantage of using a lighter target like argon is that low DM masses correspond to much lower values of v min as compared to xenon. This implies that the sensitivity forecast for a DarkSide-20k-like experiment are largely unaffected by the suppressed Gaia g(v min ) distribution at large v min .

VI. CONCLUSIONS AND OUTLOOK
The new insight into the substructure of MW's DM distribution provided by the Gaia survey forces us to move away from the simplest SHM and to re-evaluate astrophysical uncertainties in DD experiments. In this paper, we investigate the effect of Gaia Sausage, one of the most established and representative substructures, on interpreting DD data for different DM models. We demonstrate that the new Gaia velocity distribution could result in potentially large modifications of both the overall scattering rate and the recoil spectral shape. Given the limited information from existing data sets, we focus on how the Gaia velocity distribution could affect forecasting at the next generation DD experiments with the euclideanized signal method, which allows to use the full recoil spectra beyond the total recoil rate. We hope that highlighting the importance of spectral shape information could motivate the experimental collaborations to release the full three dimensional likelihood function for reconstructed DM events in case of a positive detection.
We study the sensitivity of DD experiments to different combinations of model parameters for representative DM benchmark models and its potential to distinguish different DM models given the Gaia distribution. We summarize our main findings below: • While there is still a large uncertainty in the fraction of DM in the substructure, the primary effect in reconstructing DM model parameters is due to the qualitative differences between the shapes of Gaia and SHM distributions, which are independent of the precisie value of η sub .
• For light DM with mass at or below 10 GeV, the Gaia velocity distribution leads to a significantly weakened constraint for all the models we consider. Moreover, it poses a serious challenge for identifying the DM interaction and determining its strength assuming discovery at a Darwin-like experiment with threshold E R ≈ 5 keV.
• On the contrary, for heavy DM with mass above ∼ 30 GeV, there could be a (moderate) improvement in the sensitivity of the next generation DD experiments when the mediator is heavy.
• Moreover, our results show that for positive (negative) powers of E R in the DM model, the sensitivity of a Darwin-like experiment is improved (worsened) irrespective of the DM mass and velocity distribution.
• The additional challenge in probing light DM due to the Gaia distribution could be overcome using complementary experiments with lighter targets and lower thresholds.
As mentioned in the Introduction, our work, along with several others, consists the early stages of a larger program to determine the DM phase space distribution, and assess its impact on various terrestrial experiments searching for DM. In this paper, we restrict ourselves to models with leading order SI elastic scattering for DM with mass above GeV. It will also be of interest to extend the work to models with leading order spin-dependent and/or inelastic scattering. We intend to explore how the substructures discovered using Gaia data affect the new DD experiments probing DM-electron scattering in a future publication.

Operators
Form Spin-Dependence Another important ingredient for calculating the recoil rate is the form factor, F , that encodes nuclear response functions. Ref. [65] showed that the complete basis of NREFT operators corresponds to six different types of nuclear response functions F where the subscripts indicate the NREFT operator(s). The response functions F (N,N ) i for the nuclides relevant in our analysis (see Appendix C) have been adopted from Appendix A.3 of ref. [65].

Appendix B: Euclideanized signal method
In this Appendix, we present a detailed overview of the euclideanized signal (ES) method, first introduced by refs. [41,86].
Since all DD experiments rely on counting, in some fashion, the excitations produced by DM scattering off nucleons, we can approximate the likelihood function (henceforth simply the likelihood) by a homogeneous Poisson process. For our analysis, we consider the generalized form of the likelihood given by, where d i and µ i are the number of observed and predicted events in the i th bin (the total number of bins is n b ), while the second term encodes the covariance of the background uncertainties in different bins through the matrix K ij . 11 The predicted events in each bin consist of the signal S i (θ, λ), the background B i and its associated uncertainty δB i , along with an overall normalization given by the exposure E i , Here, we have considered that the linear model is parametrized as ψ = (θ, λ, δB), such that θ and λ are the signal and nuisance parameters respectively for a given DM model. The Fisher information gives the curvature of a likelihood function around its maximum value (equivalently, the curvature around the minima of the negative log likelihood). In other words, it quantifies the maximum precision with which we can infer the model parameters from the data through the Cramér-Rao bound. Formally, the Fisher information is defined as, where E[. . .] indicates the expectation value over the observed data D(ψ). As we are working with various signal and background parameters, it's useful to represent the components of Fisher information as a block matrix, The more relevant quantity for our method is the profiled fisher information I, which can be derived using standard linear algebra rules, The individual components of the above expression can be evaluated by plugging our log-likelihood function given by eq. (B1) into eq. (B3). Simplifying further using basic algebra, we obtain where the n b × n b dimensional design matrix D is defined as, A couple of comments are in order regarding the derivation above. We have: i) assumed, following ref. [41], that the mock data is with no background perturbations, i.e δB = 0; ii) suppressed the nuisance parameters λ, and only profiled over the background uncertainties. However, these results can easily be generalized to profiling over the nuisance parameters, e.g., the substructure fraction parameter η sub . The problem of forecasting degeneracies can be interpreted as discriminating recoil spectra from two points θ 1,2 ∈ Ω M ∈ R d in a d-dimensional parameter space for a given DM model M. To this end, we can define a test statistic (TS) that can be used to reject the null hypothesis that θ 1 and θ 2 are indistinguishable at the (1 − α)% confidence level. Unfortunately, as we argue in Sec. IV, calculating the TS in eq. (26) for model comparison of N points in the parameter space can be computationally expensive when N is large. Instead, we use the ES method proposed by ref. [40], which efficiently calculates the TS by expressing it as a euclidean distance between two signals. In practice, this is achieved by transforming the parameter space into a higher-dimensional signal space such that, where ∆S = (S 1 − S 2 ) and the matrix DS is evaluated at the mean signalS = 1 2 (S 1 + S 2 ). Then, by making the embedding θ → x( θ) ∈ R n b , The TS can be written in terms of the euclidean signal (ES) x i ,  We describe configurations of the next-generation experiments used in our forecasts in this Appendix. Although there is a rich experimental program underway with different targets, we only focus on experiments that use the dual-phase (liquid-gas) time projection chamber (TPC) technology with noble elements as targets. Table IV  ton × year exposure will allow DM probes to reach the neutrino-floor for DD experiments. We model it based on the conceptual design report [6] and the latest Xenon-1T configuration. We choose the observation window in recoil energy to be  keV (cf. [100]) and divide it into 19 equally spaced bins for a primary scintillation signal (S1)-only analysis. This assumes near-perfect electron recoil (ER) background subtraction that is ensured by focusing only on the events in the nuclear recoil (NR) region. In practice, we achieve this by convolving our theory recoil spectra with the efficiency curve from the latest Xenon DM analysis (given in Fig. 1 of ref. [83]) and multiplying by a factor of 0.5. Following ref. [41], we conservatively adopt the rate for NR background components (including a 10% uncertainty) from ref. [53], rescaling them with the appropriate exposure factors. The precise choice of the background should not affect our results, since all our constraints are derived for the signal-limited region with O(100) total events.

2.
Darkside-20k (High): We follow the official proposal for a neutrino-floor LAr experiment with an integrated exposure of 200 ton×year achieved over a 10 year observation period. Our model detector has an S1-only search region with 19 linearly spaced bins in the range  keV that provides us a unique probe of O(0.1-1) TeV DM candidates. An advantage of using Ar targets is their superior pulse-shape discrimination (PSD) of the S1 signal allowing for nearly background-free detection of any DM events. Thus, we only consider a background rate of 0.1 events over the entire observation period along with a 10% systematic uncertainty. Finally, we adopt the efficiency curve for NR detection from Fig. 6 of ref. [101].
3. Darkside-20k (Low): Although there is no outlined for a low DM mass search using the Gen3 Darkside-20k setup discussed above, we make the science case for one in Sec. V D as a complementary probe to the Darwin experiment. We borrow the configuration in ref. [90] that used the Darkside-50 apparatus to perform an ionization (S2)-only analysis, albeit with a 200 ton×year exposure. There are two major differences compared to the high DM mass search: i) an S2-only analysis allows a far lower recoil energy threshold, ii) the PSD is no longer available and we need to contend with higher background rates. Thus, our fiducial analysis has a recoil energy observation window of [0. [6][7][8][9][10][11][12][13][14][15] keV and an optimistic background rate of 1 event for the entire observation period. We also use a constant acceptance of 0.43 following the discussion below Fig.1 of ref. [90]. In Fig. 12, we present constraints and forecasts in the DM coupling-mass plane for two additional models: anapole DM with a heavy mediator and magnetic dipole DM with a light mediator. They share similar features with the other DM models we discuss in detail in Sec. V A.