1 Introduction↩︎

The nature of dark matter (DM) is among the central open problems in physics. To date, evidence for DM derives from its gravitational effects on astrophysical and cosmological scales (see e.g. Refs. [1], [2] for reviews). Direct detection experiments, which search for energy depositions produced by particle DM interactions in detectors, have placed strong constraints on weak-scale DM candidates [3]–[6].

Increasing attention has therefore turned to new regions of parameter space including light sub-GeV DM, for which conventional nuclear-recoil searches can become kinematically limited [7]. In this regime electronic and collective excitations as well as other low energy processes in condensed matter targets offer promising detection channels [8]–[12]. Accessing these small energy depositions requires detectors with eV-scale or sub-eV energy thresholds, low noise and in many cases single quantum sensitivity. These capabilities are increasingly enabled by quantum sensing technologies.

Recent advances in quantum sensing have opened a new experimental window for low threshold DM searches. Representative platforms include transition edge sensors (TESs), an established class of precision calorimeters [13] that are well suited for light DM detection [10], [14]. Dedicated analysis of high resolution optical TES arrays [15] incorporating realistic detector noise and material response modeling demonstrated their potential for DM searches with thresholds down to \(\mathcal{O}(0.1)\) eV [16]. Other promising platforms include microwave kinetic inductance detectors (MKIDs) [17]–[19], among other cryogenic and non-cryogenic detector technologies. These platforms offer complementary combinations of low thresholds, energy sensitivity, scalability and low dark count rates, making them well suited to probing the small energy depositions expected from light DM interactions.

The observable particle DM interaction signals, however, depend not only on fundamental particle parameters but also on the local DM distribution and the material-dependent response of the target. Disentangling these ingredients is essential for a robust interpretation of either a potential signal or a null result. A significant challenge is that the local DM phase space distribution is not known from first principles [20]–[22]. While the Standard Halo Model (SHM) with a truncated Maxwellian velocity distribution provides a useful benchmark [23], [24], Galactic structure can include streams, substructures or even DM populations stemming from components such as a dark disk [25]–[35]. These contributions can be especially relevant for low threshold detectors, which probe small energy deposits and therefore are sensitive to low velocity regions of phase space where non-standard DM distributions can strongly affect predicted interaction rates.

Halo-independent (HI) analysis methods (see e.g. [36]–[68]), enable isolating the astrophysical dependence of direct DM detection rates by expressing them in terms of a halo function \(\widetilde{\eta}(v_{\rm min})\) where \(v_{\rm min}\) is the minimum DM speed required to produce a given recoil signal. Originally developed for nuclear recoil searches, these techniques have also been extended to DM electron scattering [69]–[71]. In the present context, the HI framework is particularly useful because detector effects can be captured by a response kernel, while the dependence on the local DM velocity distribution remains encoded in \(\widetilde{\eta}(v_{\rm min})\). This preserves a linear relation between possible measured rates and the halo function, enabling direct reconstruction and comparison across low threshold detector technologies.

In this work, we present the HI framework for sub-eV and eV energy depositions from DM scattering with electrons in detectors enabled by quantum sensing. We express the event rate in terms of the halo function and formulate a detector level linear response that incorporates finite energy resolution, threshold effects and binned data observables. As benchmarks, we analyze TES and MKID type sensors and identify the DM speed regions they probe for sub-GeV DM masses. In addition, using mock data, we demonstrate projected reconstructions of \(\widetilde{\eta}(v_{\rm min})\) for the SHM and benchmark nonstandard components, illustrating how low threshold sensors can provide complementary sensitivity to low velocity DM structures that are difficult to access with conventional recoil searches, if a putative DM signal is observed.

This paper is organized as follows. In Sec. 2, we introduce the DM electron scattering formalism relevant for quantum sensor targets. In Sec. 3, we describe the TES and MKID benchmark configurations and their material response modeling. In Sec. 4, we introduce HI formalism and define the halo function \(\widetilde{\eta}(v_{\rm min})\). In Sec. 5, we develop the detector response formalism including energy resolution and threshold effects. In Sec. 6, we define the benchmark halo scenarios, mock data event rates, binned observables, response matrix elements and show projected reconstructions for the benchmark halo functions. We conclude in Sec. 7.

2 Dark Matter Scattering↩︎

The total DM interaction rate per unit detector mass is \[R = {1\over\rho_T}{\rho_\chi\over m_\chi}\int\text{d}^3 \boldsymbol{v} f_\chi\!\left(\boldsymbol{v}\right)\Gamma\!\left(\boldsymbol{v}\right), \label{R-definition}\tag{1}\] where \(\rho_T\) is the target mass density, \(\rho_\chi\) is the local DM mass density, \(m_\chi\) is the DM particle mass, \(f_\chi\!\left(\boldsymbol{v}\right)\) is the local DM velocity distribution in Earth’s frame normalized to 1, \(\int \text{d}^3\boldsymbol{v} f_\chi(\boldsymbol{v})=1\), and \(\Gamma\!\left(\boldsymbol{v}\right)\) is the scattering rate for an incident DM particle velocity \(\boldsymbol{v}\).

In this work, we focus on DM scattering off target Standard Model (SM) electrons. Although we do not restrict ourselves to a specific underlying model, interactions of this form can readily arise in scenarios with a vector mediator of mass \(m_V\) and couplings \(g_\chi\) and \(g_e\) with the DM and electrons respectively, such as with a kinetically mixed dark photon [72], [73]. The scattering rate is \[\Gamma\!\left(\boldsymbol{v}\right)= {\pi \overline{\sigma}_e\over\mu_{\chi e}^2} \int{\text{d}^3\boldsymbol{q}\over\left(2\pi\right)^3} \left|\mathcal{F}_\text{med}\!\left(q\right)\right|^2 S\!\left(q,E_e\right), \label{Gamma32of32vecv}\tag{2}\] where \(q=\lvert\boldsymbol{q}\rvert\) denotes momentum transfer, \(\mu_{\chi e}\) is the DM-electron reduced mass, and \(S(q,E_e)\) is the material-dependent dynamic structure factor for electronic excitations. The deposited energy is determined by the non-relativistic scattering kinematics \(E_e \simeq \boldsymbol{q}\cdot\boldsymbol{v}- q^2/(2m_\chi)\), as explained in App. 8.

We define the reference DM-electron interaction cross section, evaluated at the reference momentum transfer \(q_0\), as \[\overline{\sigma}_e = {g_e^2 g_\chi^2\mu_{\chi e}^2 \over \pi\left(m_V^2+q_0^2\right)^2}. \label{refSigma}\tag{3}\] The corresponding interaction mediator form factor, that takes into account its propagator, is \[\mathcal{F}_\text{med}\!\left(q\right) = {m_V^2+q_0^2\over m_V^2+q^2}. \label{DM32form32factor}\tag{4}\] For DM-electron scattering, we adopt the conventional reference momentum \(q_0=\alpha m_e\) following the convention of Ref. [8], [69], where \(m_e\) is the electron mass and \(\alpha\) is the electromagnetic fine structure constant. For the heavy mediator limit \(\mathcal{F}_\text{med}\!\left(q\right) \simeq 1\), while for the light mediator limit \(\mathcal{F}_\text{med}\!\left(q\right) \simeq (q_0/q)^2\).

Imposing energy conservation for the transferred energy \(E\), the differential event rate is given by (see the detailed derivation in App. 9) \[{\text{d}R\over\text{d}E} = {1\over{4\pi}\rho_T\mu_{\chi e}^2} \int_{0}^{\infty} \text{d}q~ q\left|\mathcal{F}_\text{med}\!\left(q\right)\right|^2 S\!\left(q, E\right) \widetilde{\eta}\!\left(v_\text{min}\!\left(q, E\right)\right). \label{dRdE}\tag{5}\] The electronic dynamic structure factor \(S(\boldsymbol{q},E)\) depends on intrinsic properties of the target material and encodes material response to a momentum transfer \(\boldsymbol{q}\) and an energy transfer \(E\).

For interactions that couple to the electronic charge density through the longitudinal electromagnetic response of the medium, the same structure factor can be expressed in terms of the longitudinal dielectric function \(\epsilon_L(E,\boldsymbol{q})\) [70], [74] \[\label{eq:structure} S(\boldsymbol{q},E) = \frac{q^2}{2\pi\alpha} \frac{1}{1-e^{-\beta E}} \text{Im}\left[ -\frac{1}{\epsilon_L(E,\boldsymbol{q})} \right],\tag{6}\] where \(\beta=(k_{\rm B}T)^{-1}\), \(T\) is the target temperature, and \(k_{\rm B}\) is the Boltzmann constant. The quantity \(\text{Im}[-1/\epsilon_L(E,\boldsymbol{q})]\) is the energy loss function.

When a measured or computed momentum-dependent dielectric response is unavailable, the electronic response of a metallic target may be approximated in a simplified manner by that of an isotropic homogeneous electron gas. In this approximation, we use the zero-temperature Lindhard dielectric function [75]–[77] with a phenomenological damping width \(\gamma\), introduced through the replacement \(E\rightarrow E+i\gamma\), \[\epsilon_{\rm Lind}(E,q)= 1+\dfrac{3\omega_p^2}{q^2v_F^2} \Xi(E,q),\] where \[\begin{align} \Xi ={}& \dfrac{1}{2} +\dfrac{k_F}{4q} \left[ 1- \left( \dfrac{q}{2k_F} -Q_L \right)^2 \right] \ln\left[ \dfrac{ \dfrac{q}{2k_F} -Q_L+1 }{ \dfrac{q}{2k_F} -Q_L-1 } \right] \\ &+ \dfrac{k_F}{4q} \left[ 1- \left( \dfrac{q}{2k_F} +Q_L \right)^2 \right] \ln\left[ \dfrac{ \dfrac{q}{2k_F} +Q_L+1 }{ \dfrac{q}{2k_F} +Q_L-1 } \right]~, \end{align}\] where \(Q_L = (E+i\gamma)/(qv_F)\). Here, \(k_F\) and \(v_F\) are the Fermi momentum and Fermi velocity, respectively, \(\omega_p\) is the plasma frequency. Within the free-electron approximation \(v_F=k_F/m_e\), \(k_F=\sqrt{2m_eE_F}\), \(n_e=k_F^3/(3\pi^2)\), \(\omega_p=\sqrt{4\pi\alpha n_e/m_e}\), where \(E_F\) is the Fermi energy, \(n_e\) is the effective conduction electron number density.

The corresponding dynamic structure factor is obtained by substituting \(\epsilon_{\rm Lind}(E,q)\) into Eq. 6 . The finite width Lindhard model provides a transparent benchmark for the intraband response of an effective free electron medium.

3 Quantum Sensor Configurations↩︎

A broad range of quantum sensor technologies are being considered and developed for light sub-GeV dark matter searches, including calorimetric, resonant and threshold counting systems. In this work, we focus on two representative energy resolving configurations. We consider an Al target read out by TESs and TiN MKIDs. Other sensor technologies can be analyzed with a similar approach and framework once their material response, detection efficiency and measured energy response are specified.

We adopt material response models appropriate to each benchmark target. For Al, we use the tabulated dielectric response as computed with the DarkELF package [74], which combines experimental optical data with ab initio calculations. For TiN, for which a comparably established low temperature, momentum dependent energy loss function is not readily available in a suitable form for our analysis, we model the response with a free electron finite width Lindhard function described in Sec. 2.

3.1 Transition Edge Sensors↩︎

For the TES benchmark, we follow the detector configuration considered in Ref. [16] and take Al to be the active target material, corresponding to the superconducting absorber and collection volume coupled to the TES readout. The electronic response of Al is obtained from the dielectric tables of the DarkELF package [74]. These tables incorporate both conduction electron and interband contributions, retaining the material-specific response relevant to the energy and momentum transfers considered here.

The tabulated Al energy loss function we consider corresponds to a weighted superposition of Mermin oscillator contributions [78]–[80]. Each contribution is characterized by an oscillator energy, a damping width and a spectral weight. These parameters are fitted to experimental optical constants and reflection electron energy loss spectroscopy data within the CHAPIDIF framework [79], [80]. The resulting construction provides an experimentally constrained response at small momentum transfer together with a model based extension to finite momentum transfer. This captures the structure response over the energy range relevant for our analysis.

The Al energy loss function has a dominant plasmon feature near \(E \simeq 15~\text{eV}\) and lower energy interband structure, including a feature near \(E \simeq 1.5~\text{eV}\) at the upper edge of the deposited energy range that we analyze here. The plasmon frequency lies well above this range and is not relevant for the scattering interactions analyzed here. Within the analyzed window the Al response is therefore set by the structured electron-hole continuum of the data driven Mermin energy loss function, which gives the Al target richer momentum structure than a free electron model and produces distinct features of the TES response.

3.2 Microwave Kinetic Inductance Detectors↩︎

For the MKID benchmark, we consider TiN to be the active target material corresponding to both the superconducting MKID film and the absorber contributions. In contrast to Al, a suitable low temperature, momentum dependent dielectric response for TiN is not readily available in a form that can be directly incorporated into our calculation. We therefore model TiN as an effective free electron medium using the finite width Lindhard dielectric function.

A material specific TiN energy loss function analogous to the experimentally informed Al response would require dedicated optical or electron energy loss spectroscopy data together with a fit of the corresponding parameters, for example within a Mermin oscillator model [78]–[80]. In the absence of such a data set, the finite width Lindhard model provides a well defined and reproducible benchmark for the intraband electronic response.

The model is specified by the TiN mass density \(\rho_{\rm TiN}\), the effective Fermi energy \(E_F\), and the electronic damping width \(\gamma_{\rm TiN}\). We adopt the material parameters [81], [82] of \(\rho_\text{TiN} = 5.4~\)g cm\(^{-3}\), \(E_F = 3.94~\text{eV}\) and \(\gamma_\text{TiN} = 0.1 E_F\) as benchmark values. Here, \(\rho_{\rm TiN}\) determines the normalization of the event rate per unit detector mass, while \(E_F\) sets the characteristic momentum and velocity scales of the effective conduction electron population. The damping width accounts for the finite lifetime of electronic excitations and regulates the corresponding spectral features. Using these parameters, the TiN energy loss function is obtained as described in Sec. 2.

This construction captures the intraband electron-hole and plasmon response of the effective TiN conduction electron medium and provides a consistent benchmark for evaluating the MKID sensitivity. Material specific interband and crystal band contributions may introduce additional structure in the energy loss function, particularly at low deposited energies. We leave this detailed analysis for future study. The eventual inclusion of these contributions through a dedicated TiN dielectric response calculation would refine the detailed mapping between the measured spectrum and \(v_{\text{min}}\). The comparison with the experimentally informed Al response therefore also illustrates how the electronic structure of different target materials can produce complementary sensitivity to the local DM velocity distribution.

4 Halo-Independent Analysis↩︎

A HI analysis separates the dependence of the direct DM detection scattering rate on the local DM velocity distribution from the particle physics, target and detector-dependent response ingredients. The formalism was initially developed for DM-nucleus scattering [36]–[68] and was subsequently extended to DM-electron scattering [69]–[71]. The aim of HI analysis is to infer this common astrophysical function directly from measured rates, without imposing a particular functional form for the local DM distribution. This approach is complementary to the conventional halo-dependent analysis, in which a specific form of velocity distribution \(f_\chi(\mathbf{v})\) is assumed and the resulting event rates are used to derive constraints or preferred regions in the DM cross-section and mass parameter space.

For a non-relativistic DM particle of mass \(m_\chi\) depositing energy \(E_e\) through momentum transfer \(\boldsymbol{q}\), the minimum incoming speed required is obtained when \(\boldsymbol{q}\) and \(\boldsymbol{v}\) are parallel (see App. 8), with \[v_\text{min}(q, E_e) = \frac{E_e}{q} + \frac{q}{2m_\chi}. \label{eq:vmin-q}\tag{7}\] As a function of \(q\), this expression attains a unique minimum at \[q_* = \sqrt{2m_\chi E_e}, \qquad v_* = v_\text{min}(q_*, E_e) = \sqrt{\frac{2E_e}{m_\chi}}. \label{eq:qstar-vstar}\tag{8}\] Consequently, for every \(v_\text{min}> v_*\) the condition of Eq. 7 admits two momentum transfer solutions, \(q_-(v_\text{min}, E_e) < q_* < q_+(v_\text{min}, E_e)\), with details given in Eq. 28 . Both branches must be retained when changing the integration variable from \(q\) to \(v_\text{min}\), as detailed in App. 8.

For a specified DM mass and interaction model the predicted event rate can be written as a convolution of a halo function that contains the dependence on the local DM velocity distribution \(f_\chi(\mathbf{v})\) and a response function that encodes the particle physics, target material and detector response ingredients, as detailed below. We define the halo function as \[\begin{align} \widetilde{\eta}\!\left(v_\text{min}\right) &= {\rho_\chi\overline{\sigma}_e\over m_\chi} { \int\text{d}^3 \boldsymbol{v} {f_\chi\!\left(\boldsymbol{v}\right)\over v}\Theta\!\left(v - v_\text{min}\right) } \notag\\ &= {\rho_\chi\overline{\sigma}_e\over m_\chi} { \int_{v_\text{min}}^\infty \text{d}v {F\!\left(v\right)\over v} }.\label{halofuncDef} \end{align}\tag{9}\] When directional information is not retained, as in the present analysis, the astrophysical dependence can be expressed in terms of the speed \(v = |\boldsymbol{v}|\) distribution \(F(v)= v^2\int d\Omega_v f_\chi(\mathbf{v})\) that is obtained by integrating the velocity distribution over directions and does not require \(f_\chi(\mathbf{v})\) to be isotropic. Since \(F(v)\geq0\), the halo function \(\widetilde{\eta}(v_{\text{min}})\) is non-negative and non-increasing, with \(\widetilde{\eta}(\infty)=0\).

The detector frame DM velocity distribution is time dependent because of the Earth’s orbital motion around the Sun and its daily rotation. In this work, we consider only the time averaged event rate and therefore use the time averaged halo function \(\widetilde{\eta}(v_{\text{min}})\). Annual and daily modulation signals can be treated within the same HI framework but are left for future study.

The function \(\widetilde{\eta}(v_\text{min})\) that maximizes a given likelihood is piecewise constant and non-increasing, with at most \(d-1\) downward steps, where \(d\) is the number of data entries. This was first proven for a specific (unbinned) likelihood using the Karush-Kuhn-Tucker conditions [48], [62], and later clarified and generalized through convex geometry arguments [65], [67]. Equivalently, any likelihood can be maximized by a DM speed distribution of the form [67] \[\label{eq:Fdeltas} F(v)=\sum_{h=1}^{d-1} F_h~\delta(v-v_h)~,\tag{10}\] with parameters \(F_h\) and \(v_h\). This follows because any set of predicted rates can be reproduced by a distribution of this form, and the likelihood is always maximized for some set of rates. Note that, while the best-fit rates are unique, the best-fit distribution need not be.

Two sided pointwise bands in the \((v_\text{min},\widetilde{\eta})\) plane around the best fit can be defined at any chosen confidence level [62], [63], [67], although we do not compute confidence bands here. For extended likelihoods, suited to unbinned data, the best fit \(\widetilde{\eta}(v_\text{min})\) is guaranteed to be unique. For likelihoods that depend only on binned data, such as Poisson or Gaussian, it may not be the case [67]. In the latter case one defines a degeneracy band containing all functions that maximize the likelihood. When such a band is present, Wilks’ theorem does not apply and a \(\chi^2\) distribution with one degree of freedom cannot be assumed in the large sample limit [67]. Confidence levels can then be determined by Monte Carlo methods.

A key property of \(\widetilde{\eta}(v_\text{min})\) is that for a fixed DM mass and interaction model, including the mediator form factor and the reference cross-section convention, it is common to all target materials and detector technologies. Experiments operating at different detected energies and with different targets probe different windows in \(v_\text{min}\), and each detector’s response functions determine which portions of \(\widetilde{\eta}(v_\text{min})\) their data can constrain. Consequently, experiments observing the same DM signal must reconstruct mutually compatible halo functions within their overlapping sensitivity regions, up to statistical and systematic uncertainties. The HI approach therefore provides a direct means of comparing the results of several experiments without imposing a parametric model for the local DM velocity distribution.

5 Detector Response Functions↩︎

Detectors do not measure the deposited energy \(E\) directly but rather a proxy \(E'\), which may correspond to a number of phonons, photons or quasiparticles produced in the detector, which is related to \(E\) through an energy resolution function \(G(E,E')\) and a detection efficiency \(\varepsilon(E')\). The analysis is therefore carried out in terms of measurable \(E'\), and the response formalism should account for these detector effects.

We relate the recoil differential rate of Eq. 5 to the detected differential rate \[\frac{\text{d}R}{\text{d}E'}(E') = \varepsilon(E')\int_{0}^{\infty}\text{d}E~ G(E,E') \frac{\text{d}R}{\text{d}E}~. \label{eq:observed-diff-rate}\tag{11}\] For simplicity, we assume the detector efficiency can be described as \[\begin{align} \label{eq:threshold} \varepsilon(E') &= \theta \left(E' - E'_{\text{thr}}\right), \end{align}\tag{12}\] and a Gaussian energy resolution function of width \(\Delta\), \[\begin{align} G(E,E') &= \frac{1}{\sqrt{2\pi} \Delta} \exp \left[-\frac{(E-E')^2}{2\Delta^2}\right]. \end{align}\] For a TES-like and MKID-like experiments we adopt \(\Delta_{\text{TES}} = 0.036~\text{eV}\) [16] and \(\Delta_{\text{MKID}} = 0.13~\text{eV}\) [19], respectively. For simplicity, we adopt these values as deliberately broadened Gaussian standard deviation benchmarks rather than converting the reported full widths at half maximum to Gaussian standard deviations. In practice, we carry out the \(E\)-integration in Eq. 11 numerically, extending sufficiently far above the center \(E'\) of the resolution function \(G(E,E')\) for results to converge.

Our goal is to express the detected rate as a convolution in \(v_\text{min}\), \[\frac{\text{d}R}{\text{d}E'}(E') = \int_{0}^{\infty}\text{d}v_\text{min} \frac{\text{d}\mathcal{R}}{\text{d}E'}(v_\text{min}, E') \widetilde{\eta}(v_\text{min}), \label{dRdEprime1}\tag{13}\] which isolates the astrophysical quantities in \(\widetilde{\eta}(v_\text{min})\) and collects all target, particle physics and detector dependence into the kernel \(\text{d}\mathcal{R}/\text{d}E'\) [39].

We first change the integration variable in Eq. 5 from \(q\) to \(v_{\text{min}}\) using the two momentum transfer branches \(q_\pm(v_{\text{min}},E)\) described in App. 8. For fixed deposited energy \(E\), the map from \(q\) to \(v_{\text{min}}\) is two to one for \(v_{\text{min}}>v_*(E)\) (see Eq. 8 ). Splitting the momentum integral at \(q_*=\sqrt{2m_\chi E}\) and including the Jacobian of each branch gives

The absolute values of the Jacobians account for the opposite orientations of the two branches with \(\partial q_-/\partial v_{\text{min}}<0\), whereas \(\partial q_+/\partial v_{\text{min}}>0\).

At fixed \(v_{\text{min}}\), the deposited energy is bounded by \[E\leq E_{\text{max}}(v_{\text{min}})=\tfrac{1}{2}m_\chi v_{\text{min}}^2, \label{eq:Emax}\tag{15}\] which follows from requiring the momentum-transfer solutions to be real, \(m_\chi^2 v_{\text{min}}^2 - 2m_\chi E \geq 0\). Thus \(E_{\text{max}}(v_{\text{min}})\) is the largest energy a DM particle of minimum speed \(v_{\text{min}}\) can deposit considering optimal momentum transfer. Material excitation thresholds are already encoded in the dynamic structure factor \(S(q,E)\).

Inserting the recoil rate \(\text{d}R/\text{d}E\) of Eq. 14 into Eq. 11 and exchanging the order of the \(E\) and \(v_\text{min}\) integrations, the detected rate becomes

where the deposited energy integration extends up to the kinematic limit \(E_{\text{max}}(v_\text{min})\) of Eq. 15 , and the exchange of integration order maps the lower limit \(v_*(E)\) of Eq. 14 onto this upper limit on \(E\). Using Eqs. 13 and 16 , one can now identify the differential response function \[\begin{align} {\text{d}\mathcal{R}\over\text{d}E^\prime}\!\left(v_\text{min}, E^\prime\right) &= \sum_{\lambda=\pm}{\text{d}\mathcal{R}_\lambda\over\text{d}E^\prime}\!\left(v_\text{min}, E^\prime\right), \end{align}\] with the per branch kernels \[\begin{align} {\text{d}\mathcal{R}_{\pm}\over\text{d}E^\prime}\!\left(v_\text{min}, E^\prime\right) =&~ \frac{1}{4\pi\rho_T\mu_{\chi e}^2} \varepsilon\!\left(E^\prime\right) \int_{0}^{E_{\text{max}}(v_\text{min})}\text{d}E G\!\left(E, E^\prime\right) \nonumber\\ & \times \left|\partial q_\pm\over\partial v_\text{min}\right| q_\pm\!\left(v_\text{min}\right) \left|\mathcal{F}_\text{med}\!\left(q_\pm\!\left(v_\text{min}\right)\right)\right|^2 \nonumber\\ & \times S\!\left(q_\pm\!\left(v_\text{min}\right), E\right). \end{align}\]

In Figs. 1 and 2 we show the differential response function \(\text{d}\mathcal{R}/\text{d}E^\prime\) as a function of \(v_\text{min}\) for the TES and MKID configurations for DM electron scattering in the heavy mediator limit and light mediator limit, respectively. The response functions act as window functions in \(v_\text{min}\). Measurements at a given observed energy \(E^\prime\) constrain the halo function only over the range of \(v_\text{min}\) for which the corresponding response are different from zero, shown more clearly in the bottom panels of both figures, with the function normalized to a maximum of approximately 1. The lower edge of this range is not determined solely by the detector threshold. The threshold enters through the efficiency \(\epsilon(E^\prime)\) in Eq. 12 , whereas the \(v_\text{min}\) dependence follows from the kinematic minimum \(v_*\) evaluated over the true deposited energies \(E\) selected by the resolution function \(G(E,E^\prime)\) and weighted by the material response \(S(q,E)\). The finite energy resolution therefore smooths and broadens the mapping between \(E^\prime\) and \(v_\text{min}\).

For both TES and MKID detector configurations we consider, increasing \(E'\) generally shifts the response toward larger \(v_\text{min}\), while increasing the DM mass shifts the probed range toward smaller \(v_\text{min}\), as can be approximately noted from the scaling \(v_*\propto m_\chi^{-1/2}\). The detailed features are target dependent. The Al response used for TES contains several interband and collective excitation features, which produce the multiple peaks visible in Fig. 1, particularly for \(m_\chi=1~\text{GeV}\) scenario. The effective free electron TiN response used for MKID is comparatively smoother and selects a different, although partially overlapping, range of \(v_\text{min}\). The bottom panels of Fig. 1 and Fig. 2, in which each response is normalized to its maximum, make these distinct velocity windows and their complementarity more apparent.

For the light mediator case since \(\mathcal{F}_\text{med}(q) \propto q^{-2}\) the resulting strong weighting toward small momentum transfer suppresses the large \(q\) branch relative to the low \(q\) branch and thus redistributes the response in \(v_\text{min}\). Consequently, both the normalization and the shape of the \(v_{\rm min}\) space window depend on the assumed mediator behavior. The HI reconstruction is therefore model independent with respect to the DM velocity distribution, but depends on the chosen DM mass and interaction model.

6 Benchmark Halo Reconstruction↩︎

In this section we apply the method developed above to mock data, using the response functions obtained in Sec. 5. We generate event rates for several benchmark local DM velocity distribution models, adopting a realistic local DM density \(\rho_\text{DM} \simeq 0.4~\text{GeV}\) \(^{-3}\) (see e.g. Ref. [83]) together with a deliberately large reference DM-electron cross section of \(\sigma_e = 10^{-30}~\text{cm}^2\). The HI method is most informative when a sufficient number of signal events is observed to reconstruct an energy spectrum. Hence, to illustrate this, the large cross section is chosen so that the predicted event counts are sufficiently high to expose the performance of the reconstruction method itself. We note that this cross section choice is illustrative. For our analysis we treat the resulting mock event rates as “observed” data and reconstruct the input halo function used to generate them, assuming a heavy mediator throughout.

From the analysis we obtain piecewise constant best-fit halo functions as expected. We find that TES-like and MKID-like example experiments yield complementary results.

6.1 Dark Matter Velocity Distribution Models↩︎

We analyze the reconstruction of three benchmark local DM velocity distributions. The first is the Standard Halo Model (SHM), which serves as our baseline. The second is the SHM with an added co-rotating dark disk (DD). The third is the SHM with an added Earth-bound (EB) population. The two non-standard benchmarks beyond pure SHM introduce low \(v_\text{min}\) structure of the kind that low threshold quantum sensors are particularly well equipped to probe. We define the three models and their halo functions \(\widetilde{\eta}(v_\text{min})\) below and quote the closed form expressions in App. 10. The corresponding curves and their reconstructed best fits are shown in Figs. 6, 7 and 8.

For the SHM the local distribution is an isotropic Maxwell-Boltzmann in the Galactic frame, truncated at the Galactic escape speed and boosted to the detector frame by Earth’s motion. In terms of the DM velocity \(\boldsymbol{v}\) in the lab frame velocity, with the lab at rest with respect to Earth, \[\begin{align} f_\chi(\boldsymbol{v}) = \frac{1}{N_0}~ e^{-(\boldsymbol{v} + \boldsymbol{v}_e)^2/v_0^2}~ \Theta\!\left(v_\text{esc} - |\boldsymbol{v} + \boldsymbol{v}_e|\right), \label{SHM95dif} \end{align}\tag{17}\] where \(N_0\) normalizes the distribution to unity, \(v_0\) is the local circular speed with velocity dispersion \(\sigma_0 = v_0/\sqrt{2}\), \(\boldsymbol{v}_e\) is Earth’s mean velocity relative to the Galaxy and \(v_\text{esc}\) is the Galactic escape speed at the solar location. We adopt the conventional values \(v_0 = 238~\)km s\(^{-1}\), \(v_e = 250~\)km s\(^{-1}\) and \(v_\text{esc} = 544~\)km s\(^{-1}\) [84]. With these values \(\widetilde{\eta}_\text{SHM}(v_\text{min})\), given in closed form in App. 10, varies slowly at low \(v_{\rm min}\), changes analytic form at \(v_\text{esc}-v_e \simeq 294~\)km s\(^{-1}\), and vanishes when \(v_\text{min}\geq v_\text{esc}+v_e \simeq 794~\)km s\(^{-1}\) as displayed in Figs. 6, 7 and 8.

As a first benchmark beyond SHM we consider an additional contribution from a DD. Such a component can arise as a relic of the Galaxy’s mergers history or in models where a subdominant DM fraction undergoes dissipative self-interactions and settles into a rotationally supported structure, which can be approximately co-planar with the baryonic disk and lagging behind it in rotation [25]–[29]. In the Solar neighborhood this results in a small dispersion and a modest bulk velocity relative to Earth and an excess in \(\widetilde{\eta}(v_\text{min})\) at low \(v_\text{min}\), where the SHM is nearly flat.

We adopt the DD model of Ref. [85], a truncated Maxwellian form with \(v_0^\text{DD} = 70~\)km s\(^{-1}\), \(v_e^\text{DD} = 100~\)km s\(^{-1}\) (which is Earth’s speed with respect to the DD rest frame) and \(v_\text{esc}^\text{DD} = 694~\)km s\(^{-1}\) (which is the escape speed with respect to the DD). Here, both dispersion and bulk speed are well below their SHM counterparts, so that \(\widetilde{\eta}_\text{DD}(v_\text{min})\) is concentrated at low \(v_\text{min}\). The total halo function is the density weighted combination \[\begin{align} \widetilde{\eta}_\text{SHM+DD}(v_\text{min}) = (1 - f_\text{DD})~\widetilde{\eta}_\text{SHM}(v_\text{min}) + f_\text{DD}~\widetilde{\eta}_\text{DD}(v_\text{min}), \label{etaSHMpDD} \end{align}\tag{18}\] where \(f_\text{DD}\) is the DD fraction of the local density \(\rho_\text{DM} = 0.4~\text{GeV}\) \(^{-3}\) and \(\widetilde{\eta}_\text{DD}\) follows from Eq. 9 with the DD distribution. Although DD fractions as large as \(\sim25\%\) were previously considered [27], more recent analyses following Gaia survey observations disfavor values above \(\sim5\%\) [30]–[32]. We display both \(f_\text{DD} = 5\%\) and \(25\%\) cases as illustrative reconstructions.

As a second, more extreme benchmark we add a gravitationally Earth-bound population, which can build up as DM particles lose energy in repeated scatterings in the Earth or Sun and become captured in bound orbits [33]–[35]. Earth capture proceeds through DM scattering on SM constituents such as baryons, thus the capture rate and the resulting bound density are model-dependent. We do not link this capture to the DM-electron scattering interaction that we analyze as producing signals in experiments. Instead we treat the bound DM population phenomenologically remaining agnostic regarding its specific model-dependent particle interaction origin, fixing its density and velocity parameters directly and thus isolating its effect on the reconstruction. The assumed bound DM density can be many orders of magnitude larger than the ambient local DM density. Its very low lab frame velocities produce a steep enhancement of \(\widetilde{\eta}(v_\text{min})\) as \(v_\text{min}\to 0\), which for sub-GeV DM can be accessed by sub-eV threshold quantum sensors.

Following Ref. [85] we model this scenario as a truncated Maxwellian of the form of Eq. 17 with values \(\rho_\text{EB} = 10^{14}~\text{GeV}\) \(^{-3}\), \(v_0^\text{EB} = \sqrt{2k_BT/m_\chi}\), \(v_e^\text{EB} = 0\), \(v_\text{esc}^\text{EB} = 11.2~\)km s\(^{-1}\) [33]–[35]. Here, \(T \simeq 300~\text{K}\) is the ambient terrestrial temperature, so that \(v_0^\text{EB}\) is the thermal speed in equipartition with terrestrial matter, so that \(v_0^\text{EB} = 2.16~\)km s\(^{-1}\) for \(m_\chi = 1~\text{GeV}\), and \(v_e^\text{EB} = 0\) encodes the co-moving character of a bound component, and \(v_\text{esc}^\text{EB}\) is Earth’s surface escape speed. The total halo function is additive, \[\begin{align} \widetilde{\eta}_\text{SHM+EB}(v_\text{min}) = \widetilde{\eta}_\text{SHM}(v_\text{min}) + \widetilde{\eta}_\text{EB}(v_\text{min}), \label{etaSHMpEB} \end{align}\tag{19}\] each term carrying its own density prefactor in Eq. 9 with \(\rho_\text{DM}\) for the SHM and \(\rho_\text{EB}\) for the EB. Since \(\rho_\text{EB}/\rho_\text{DM} \simeq 10^{14}\), the EB term dominates \(\widetilde{\eta}_\text{SHM+EB}\) throughout the accessible region \(v_\text{min}\lesssim v_\text{esc}^\text{EB}\).

6.2 Mock Data Event Rates↩︎

Although the detector configurations considered here are sensitive to individual events, reconstructing the halo function benefits from spectral information. We therefore analyze the distribution of assumed signal events in bins of observed energy. For concreteness, for both detector configurations we analyze we bin the observed energy \(E'\) into five bins spanning from each detector’s threshold up to \(1.1~\mathrm{eV}\). The upper value is chosen to concentrate on the sub-eV to eV window enabled by low threshold quantum sensors, which is where these detectors provide reach beyond conventional direct DM detection experiments and which maps onto the low \(v_\text{min}\) region we aim to reconstruct. For the TES configuration, with a threshold of \(0.1~\mathrm{eV}\), these are five uniform bins of width \(0.2~\mathrm{eV}\) giving \([0.1,0.3]\), \([0.3,0.5]\), \([0.5,0.7]\), \([0.7,0.9]\) and \([0.9,1.1]~\mathrm{eV}\). For the MKID configuration, whose threshold is \(0.2~\mathrm{eV}\), we keep the same upper four bins and use a narrower lowest bin \([0.2,0.3]~\mathrm{eV}\), giving \([0.2,0.3]\), \([0.3,0.5]\), \([0.5,0.7]\), \([0.7,0.9]\) and \([0.9,1.1]~\mathrm{eV}\).

In the \(i\)-th observed energy bin \([E_i^\prime, E_{i+1}^\prime]\), the event rate is \[\begin{align} R_i = R_{[E_i^\prime, E_{i+1}^\prime]} = \int_{E_i^\prime}^{E_{i+1}^\prime}\text{d}E^\prime {\text{d}R\over\text{d}E^\prime}\!\left(E^\prime\right), \end{align}\] and the expected predicted number of events in that bin is \[\begin{align} \braket{N_\text{event}}_i = M_T T R_i, \label{eventNumber} \end{align}\tag{20}\] where \(M_T\) is the detector mass and \(T\) the total running time. For the projected exposures, we adopt \[\begin{align} M_T T = \begin{cases} 8.2~\mathrm{\mu g\cdot month} & (\text{TES}~\cite{Chen:2025cvl}),\\[2pt] 10^7~\mathrm{pixel\cdot year} = 4.2~\mathrm{m g\cdot year} & (\text{MKID}~\cite{Gao:2024irf}), \end{cases} \label{exposures} \end{align}\tag{21}\] where a target mass of \(0.42~\text{ng}\) per pixel is considered for the MKID configuration [86]. The computed expected event numbers for all three benchmark models at these exposures are shown in Figs. 3 and 8.

6.3 Reconstruction of Speed Distributions↩︎

The event rate in the \(i\)-th observed energy bin can be stated as a convolution of the halo function with the response integrated over that bin, \[\begin{align} R_i = \int \text{d}v_\text{min} \mathcal{R}_{\left[E_i^\prime, E_{i+1}^\prime\right]}\!\!\left(v_\text{min}\right) \widetilde{\eta}\!\left(v_\text{min}\right), \end{align}\] where the integrated response function is \[\begin{align} \mathcal{R}_{\left[E_i^\prime, E_{i+1}^\prime\right]}\!\!\left(v_\text{min}\right) = \int_{E_i^\prime}^{E_{i+1}^\prime}\text{d}E^\prime {\text{d}\mathcal{R}\over\text{d}E^\prime}\!\left(v_\text{min}, E^\prime\right). \label{binCurlyR} \end{align}\tag{22}\] This function is shown in Fig. 4 for the TES and MKID benchmarks over several \(E^\prime\) bins. A given bin constrains \(\widetilde{\eta}(v_\text{min})\) only over the range of \(v_\text{min}\) where its \(\mathcal{R}_{\left[E_i^\prime, E_{i+1}^\prime\right]}\) is appreciably nonzero. Thus, the integrated responses effectively act as window functions in \(v_\text{min}\).

Following the procedure of Ref. [47], we discretize the halo function on the \(v_\text{min}\) interval of interest, namely in the range over which the response functions are appreciably nonzero as shown in Fig. 4. We divide it into \(N_\text{int}\) equal adjacent intervals \(j = 1,\dots,N_\text{int}\) on each of which \(\widetilde{\eta}\) takes a unique value \(\widetilde{\eta}_j\), \[\begin{align} \widetilde{\eta}_{\rm ansatz}(v_\text{min}) = \sum_{j=1}^{N_\text{int}} \widetilde{\eta}_j \left[\Theta\!\left(v_\text{min}- v_\text{min}^j\right) - \Theta\!\left(v_\text{min}- v_\text{min}^{j+1}\right)\right], \label{etaStepForm} \end{align}\tag{23}\] as illustrated in Fig. 5. The \(\widetilde{\eta}_j\) are the parameters fitted by the constrained minimization described below. In the best-fit, multiple adjacent intervals take the same value. Hence, the reconstructed \(\widetilde{\eta}(v_\text{min})\) is piecewise constant with at most \((d-1)\) downward steps, where \(d\) is the number of data bins. The interval boundaries serve to locate these steps, since each step must fall on one of them. We therefore use a grid much finer than the expected number of steps and verify convergence to a particular best-fit function. Namely starting from a number of intervals larger than the number of data bins and doubling it, the best fit shape ceases to change once the grid is sufficiently fine to resolve the step locations, after which further refinement leaves them fixed.

Inserting the ansatz Eq. 23 into the rate, the predicted differential rate is \[\begin{align} {\text{d}R\over\text{d}E^\prime}\!\left(E^\prime\right) = \sum_{j=1}^{N_\text{int}} \widetilde{\eta}_j \int_{v_\text{min}^j}^{v_\text{min}^{j+1}}\text{d}v_\text{min} {\text{d}\mathcal{R}\over\text{d}E^\prime}\!\left(v_\text{min}, E^\prime\right). \end{align}\] Its integral over bin \(i\) is \[\begin{align} R_i =&~ \sum_{j=1}^{N_\text{int}} \widetilde{\eta}_j \int_{E_i^\prime}^{E_{i+1}^\prime}\text{d}E^\prime \int_{v_\text{min}^j}^{v_\text{min}^{j+1}}\text{d}v_\text{min} {\text{d}\mathcal{R}\over\text{d}E^\prime}\!\left(v_\text{min}, E^\prime\right) \nonumber\\ =&~ \sum_{j=1}^{N_\text{int}} \mathcal{R}_{ij}\widetilde{\eta}_j, \end{align}\] which defines the response matrix elements \[\begin{align} \mathcal{R}_{ij} =&~ \int_{E_i^\prime}^{E_{i+1}^\prime}\text{d}E^\prime \int_{v_\text{min}^j}^{v_\text{min}^{j+1}}\text{d}v_\text{min} {\text{d}\mathcal{R}\over\text{d}E^\prime}\!\left(v_\text{min}, E^\prime\right) \nonumber\\ =&~ \int_{v_\text{min}^j}^{v_\text{min}^{j+1}}\text{d}v_\text{min} \mathcal{R}_{\left[E_i^\prime, E_{i+1}^\prime\right]}\!\!\left(v_\text{min}\right), \label{responseMatrix} \end{align}\tag{24}\] with the second form following from Eq. 22 .

As mock observed data \(\mathcal{O}_i\) we take the per bin event numbers \(\mathcal{O}_i = \braket{N_\text{event}}_i\) of Eq. 20 , computed for each of the three benchmark DM velocity distributions. We obtain the best fit \(\widetilde{\eta}(v_\text{min})\) by minimizing the Neyman \(\chi^2\) statistic, \[\begin{align} \chi^2_\text{N} = \sum_i \frac{(\mu_i - \mathcal{O}_i)^2}{\mathcal{O}_i}, \label{chiSq} \end{align}\tag{25}\] in which the variance of each bin is approximated by its observed Poisson count \(\mathcal{O}_i\), so that the statistical uncertainty is \(\sigma_{\text{stat},i} \simeq \sqrt{\mathcal{O}_i}\). This Gaussian approximation to the Poisson likelihood is accurate for bins with about ten events or more. It is one convenient choice and not a strict requirement of the HI method, and other statistics such as the Poisson likelihood \(\chi^2\) of Ref. [87] could be used instead. In this simplified analysis we neglect backgrounds, such as a possible low energy excess [88], as well as systematic uncertainties.

The predicted event number per bin is \[\begin{align} \mu_i = M_T T \sum_{j=1}^{N_\text{int}} \mathcal{R}_{ij} \widetilde{\eta}_j, \label{mudef} \end{align}\tag{26}\] and the minimization procedure determines \(\widetilde{\eta}_j\), yielding an approximation to the input \(\widetilde{\eta}\) function. During the minimization we impose two physical constraints, \[\begin{align} \widetilde{\eta}_j \geq 0, \qquad \widetilde{\eta}_j \geq \widetilde{\eta}_{j+1} \quad (j = 1,\dots,N_\text{int}-1), \end{align}\] enforcing the positivity of the DM speed distribution and the non-increasing character of any physical halo integral.

The statistical Poisson uncertainties on the input counts \(\mathcal{O}_i\) are shown in Fig. 3. Propagating them through the constrained minimization requires Monte Carlo methods and is deferred to future work.

We note that since our likelihood is based on binned event counts, as explained in Sec. 4, the best fit halo function need not be unique [67]. The minimization therefore returns a piecewise constant best-fit solution. Other halo functions may yield the same minimum value of the test statistic and identical predicted bin counts. We do not determine the corresponding degeneracy or confidence bands here.

6.4 Reconstructed Halo Functions↩︎

In Figs. 6 7 and 8 we show the original and best-fit \(\widetilde{\eta}(v_\text{min})\) functions, for TES and MKID configurations, for \(m_\chi = 10\) MeV and 1 GeV.

In the analysis, for each mass and detector configuration the \(v_\text{min}\) axis is divided into \(N_{\rm int}\) intervals covering the kinematically accessible window. Concretely, for the TES (Al) and \(m_\chi= 10\) MeV, we divide the window \(\left[42.5~\text{km}~\text{s}^{-1}, 150~\text{km}~\text{s}^{-1}\right]\) into \(135\) intervals, and for \(m_\chi= 1\) GeV, we divide \(\left[4.2~\text{km}~\text{s}^{-1}, 78~\text{km}~\text{s}^{-1}\right]\) into \(92\) intervals (see in Fig. 4 where the response functions are non-zero). For the MKID (TiN), we divide \(\left[56.1~\text{km}~\text{s}^{-1}, 168~\text{km}~\text{s}^{-1}\right]\) into \(140\) intervals when \(10\;\mathrm{MeV}\) and \(\left[5.0~\text{km}~\text{s}^{-1}, 134~\text{km}~\text{s}^{-1}\right]\) into \(162\) intervals when \(1\;\mathrm{GeV}\). In every case the intervals are equally spaced, with width \(0.8~\text{km}~\text{s}^{-1}\). Since the spacing is linear while the features that distinguish the benchmark models lie at low \(v_\text{min}\), a small fraction of the linear range but a broad region on the logarithmic \(v_\text{min}\) axis of Figs. 6, 7 and 8, a large \(N_\text{int}\) is required to position sufficient interval boundaries at low velocities. The height of the horizontal bars indicate constant pieces of the best-fit, which are at most five, and their projection into the horizontal axis the \(v_\text{min}\) interval they correspond to.

In Fig. 6, for \(m_\chi = 10~\text{MeV}\), the \(v_\text{min}\) ranges to which both detector configurations are sensitive correspond to the \(\widetilde{\eta}_{\rm SHM}\) slowly varying region, below where \(\widetilde{\eta}_\text{SHM}\) falls near \(v_\text{min}\simeq 294~\)km s\(^{-1}\). Their reconstructions closely trace the input. TES and MKID cover offset although partially overlapping \(v_\text{min}\) windows, determined by their different material responses (i.e. Al versus TiN) and resolutions. In the overlapping range the two best fits agree with each other and with the SHM, which serves as an internal consistency check of the method since a disagreement there could point to detector systematic effects or an incorrect response model rather than to a real feature of \(\widetilde{\eta}\). For \(m_\chi = 1~\text{GeV}\) the kinematic minimum \(v_{\ast}\) is smaller by a factor \(\mathcal{O}(10)\), and both detectors reach \(\mathcal{O}(10)~\)km s\(^{-1}\). There \(\widetilde{\eta}_\text{SHM}\) is essentially constant, thus the reconstruction primarily fixes the normalization factor \(\rho_\text{DM} \sigma_e/m_\chi\) rather than the spectral shape and the best fits remain mutually consistent and stable even as the response kernel becomes more sharply peaked at low \(v_\text{min}\) as can be seen in Fig. 4.

The SHM reconstruction analysis thus establishes two baseline properties. First, the reconstruction recovers the slowly varying portion of the input across the accessible \(v_\text{min}\) range at both considered DM masses. Second, TES and MKID detector configurations already display complementarity through their offset but overlapping windows, and hence that combined measurements between experiments can both broaden the reconstructed halo range and provide consistency checks in the overlapping range. By construction the SHM has no significant velocity-dependent structure in this window and the discriminating power of the method could highlight possible components additional to the SHM.

In Fig. 7 we display the SHM+DD reconstructions for the two masses and two contributing DM disk density fractions. In each panel the solid curve is the full SHM+DD prediction of Eq. 18 , the dashed curve depicts pure SHM and the dash-dotted curve depicts the SHM scaled by \((1-f_\text{DD})\) corresponding to the halo function at fixed total density with the disk component removed. The gap between the dash-dotted and solid curves at low \(v_\text{min}\) is the isolated disk contribution \(f_\text{DD} \widetilde{\eta}_\text{DD}\), hence the DD signature is directly readable from the bin distribution.

For \(f_\text{DD} = 5\%\) case the departure from the pure SHM is found to be small across the accessible window range, and the best fit bins follow the SHM+DD curve while remaining close to the SHM behavior. Thus, discriminating quantitatively between the two possibilities would require analysis of the uncertainty bands that is beyond our scope. For the more extreme case of \(f_\text{DD} = 25\%\) the excess is clearly more pronounced. At \(m_\chi = 10~\text{MeV}\) the accessible window lies near the peak of \(\widetilde{\eta}_\text{DD}\), and the low \(v_\text{min}\) best-fit TES bins rise above the SHM curve, accommodating DD contributions, onto the SHM+DD curve. For \(m_\chi = 1~\text{GeV}\) the window shifts further into the range dominated by disk contributions, at lower \(v_\text{min}\). The MKID configuration, with its higher energy threshold and different TiN response, covers a higher \(v_\text{min}\) window and anchors the transition where SHM and SHM+DD curves converge. Together, the detector configurations constrain both the amplitude of the excess beyond SHM and the speed range where it drops, determined by \(v_0^\text{DD}\) and \(v_e^\text{DD}\). The framework therefore could enable, subject to detailed analysis of uncertainties, potentially resolving low \(v_\text{min}\) features once their amplitudes are comparable to or larger than the SHM plateau. The detector complementarity is shown to play a prominent role especially when the features in halo function are less pronounced or the accessible window ranges are narrow.

In Fig. 8 we display the SHM+EB reconstruction at \(m_\chi = 1~\text{GeV}\) on a logarithmic scale for \(\tilde{\eta}\). We restrict to \(1~\text{GeV}\) since \(v_{\ast}\) must fall below \(v_\text{esc}^\text{EB} \simeq 11.2\) km s\(^{-1}\) for the EB component to noticeably contribute, and at \(10~\text{MeV}\) this lies below the detector thresholds leaving no accessible bins sensitive to EB components. Reaching lighter masses requires detector configurations with lower thresholds. The model curve rises by about a decade in \(\widetilde{\eta}\) per decade in \(v_\text{min}\) as \(v_\text{min}\to 0\), reflecting the dramatic potential EB overdensity beyond SHM in the region below Earth’s escape speed. Both detector configurations we consider are capable of recovering this steep rise within their window ranges without any assumed functional form of the halo model. This highlights regimes where the HI approach yields information unavailable to conventional direct DM detection parametric halo fits. The right panel shows the corresponding per bin event numbers that peak in the lowest accessible \(E'\) bins, the observational counterpart of the low \(v_\text{min}\) enhancement mapped through \(\mathcal{R}_{ij}\). We stress that the very large counts follow directly from the assumed \(\rho_\text{EB} \simeq 10^{14}\) GeV cm\(^{-3}\) overdensity combined with \(\sigma_e = 10^{-30}~\text{cm}^2\), which is not realistic and depicted for illustration.

Across the three considered halo function benchmarks that reflect scenarios respectively with a flat plateau, a modest low \(v_\text{min}\) excess and a steep decade per decade rise, the HI reconstruction is seen to appropriately track the input \(\widetilde{\eta}\) functions and shows consistent complementarity of detector configurations. In each case the analysis output is a best-fit piecewise constant \(\widetilde{\eta}\) function, obtained without assuming a specific halo ansatz, which carries two key physical information components, related to halo distribution normalization and spectral shape. The normalization is associated with prefactor \(\rho_\chi\sigma_e/m_\chi\) of Eq. 9 . For fixed \(m_\chi\) and \(\rho_\chi\), the plateau level of the best fit \(\widetilde{\eta}\) function constrains \(\sigma_e\) and hence the coupling product \(g_e g_\chi\) that enters the DM-electron amplitude, even if the velocity distribution departs from the SHM. The spectral shape across bins indicates the underlying model scenario that is nearly flat for the SHM as shown in Fig. 6, has a modest low \(v_\text{min}\) excess for a DD as shown in Fig. 7 and has a steep rise toward \(v_\text{min}\to 0\) for an Earth-bound population as shown in Fig. 8. The best fits already clearly illustrate different low \(v_\text{min}\) trends across the three scenarios, however quantitative model discrimination requires detailed analysis of uncertainty bands on the \(\widetilde{\eta}_j\) that we defer to future work.

The HI reconstruction also carries an incorporated restriction for physical results. Any \(f_\chi(v)\) yields a non-increasing \(\widetilde{\eta}(v_\text{min})\). Hence, the monotonicity constraint described above is automatically satisfied by a genuine DM signal. This restriction requires no additional astrophysical input and is present in conventional direct DM detection parametric halo analysis fits. Moreover, because \(\widetilde{\eta}(v_\text{min})\) is detector independent a reconstruction from one experiment can be propagated to predict spectrum in another experiment. Inserting results from TES-derived \(\widetilde{\eta}\) analysis into Eq. 26 with the MKID response matrix enables predicting the MKID experiment counts, and vice versa. Thus, the HI method enables a direct cross experiment consistency test as the halo function reconstructed from one experiment can be used to predict the spectrum in another, and agreement with the observed spectrum can be assessed without assuming a specific halo model.

We further comment on the origin of the detector complementarity. The Al target used for the TES benchmark is modeled with a data driven Mermin energy loss function that includes material specific interband and collective structure. By contrast the TiN target used for the MKID benchmark is described by a smoother finite width Lindhard response containing only the intraband free electron contribution. These differences lead to distinct response functions and sensitivity to different \(v_{\text{min}}\) ranges. Neither the Al nor the TiN plasmon enters the sub-eV to eV window analyzed here as the Al plasmon resides near \(\sim 15~\text{eV}\) and the free-electron TiN plasmon is several eV above the window. Hence, in both cases the response is governed by the low energy electron-hole continuum. Through \(v_\text{min}\), the richer momentum structure of the Al energy loss function produces the multiply peaked TES response of Fig. 1. For TiN the response is comparatively smoother. Hence, the two detector configurations probe structurally different parts of \(\widetilde{\eta}(v_\text{min})\).

Since binned data may admit degenerate best fit halo functions, a statistically robust cross-experiment prediction requires propagating the corresponding degeneracy or confidence set, or performing a joint fit. In the present analysis we illustrate the mapping using only the best-fit solutions.

6.5 Other Detector Applications↩︎

The HI formalism described here applies to any sub-eV sensor once its response function is computed. Notably, this is not restricted to energy resolving detectors. A threshold (i.e. counting) detector that does not resolve the deposited energy spectrum still contributes a well defined response function that is determined by its threshold and target material. If different thresholds can be achieved they could be combined into binned data.

One such example is superconducting nanowire single photon detectors (SNSPDs), reaching sub-eV thresholds with low dark count rates while their superconducting films (e.g. based on NbN or WSi) provide additional target materials. The QROCODILE concept, based on SNSPD readout, reaches very low thresholds and can extend \(v_\text{min}\) coverage [89]. Energy resolving platforms such as SQUID-read magnetic microcalorimeters and bolometers provide a complementary spectrometric handle. Hence, a multi-target, multi-threshold program combining these channels is a natural extension of the low \(v_\text{min}\) HI approach.

7 Conclusions↩︎

We have presented a HI framework for interpreting sub-GeV DM signals in quantum sensors with sub-eV thresholds. By separating the common halo function \(\widetilde{\eta}(v_{\text{min}})\) from the particle physics, material and detector response, the method enables translating measured experimental spectra into information about the local DM velocity distribution without assuming a specific halo model.

Using TES (Al) and MKID (TiN) benchmarks, we have shown that different sensor materials and thresholds probe complementary, partially overlapping regions of \(v_{\text{min}}\). Their combination can therefore extend the velocity range accessible to reconstruction and provide a direct consistency test. A genuine DM signal must yield a common halo function across detectors. The benchmark SHM, dark disk and Earth-bound scenarios further illustrate how non-standard low velocity populations could leave distinguishable features in the reconstructed distribution. The detector response functions used here incorporate realistic material and detector inputs. The reconstructions themselves are intended as demonstrations of the method and employ deliberately high signal statistics. Establishing realistic discovery and discrimination sensitivity reach will require detailed confidence bands and a complete treatment of statistical and systematic uncertainties, which we leave for future work.

More broadly, this framework turns low threshold quantum sensors from instruments that only constrain interaction strengths into potential probes of the local DM environment. It can be extended to other energy resolving or threshold counting platforms once their response functions are known. This motivates a multi-material and multi-platform program capable of mapping DM velocities in a regime that is difficult to access with conventional direct DM detection experiments.

Acknowledgments↩︎

G.B.G.is partially supported by the US Department of Energy under Award Number DE-SC0009937. V.T.and M.C.were supported by the World Premier International Research Center Initiative (WPI), MEXT, Japan. V.T.acknowledges support from JSPS KAKENHI grant No. 23K13109.

8 General Scattering Formalism↩︎

Consider a non-relativistic DM particle of mass \(m_\chi\) with incoming momentum \(\boldsymbol{p}=m_\chi\boldsymbol{v}\) and outgoing momentum \(\boldsymbol{p}^\prime=\boldsymbol{p}-\boldsymbol{q}\), where \(\boldsymbol{q}\) is the momentum transferred to the target.

Energy conservation implies that the deposited energy \(E_e\) is \[\begin{align} E_e =E_i-E_f =\frac{\boldsymbol{p}^2}{2m_\chi} -\frac{(\boldsymbol{p}-\boldsymbol{q})^2}{2m_\chi} =\frac{2\boldsymbol{p}\cdot\boldsymbol{q}-q^2}{2m_\chi}, \end{align}\] where \(E_i=\boldsymbol{p}^2/(2m_\chi)\) and \(E_f=\boldsymbol{p}^{\prime 2}/(2m_\chi)\) are the initial and final DM kinetic energies. Thus, \[\begin{align} E_e &=\boldsymbol{q}\cdot\boldsymbol{v}-\frac{q^2}{2m_\chi} =qv\cos\theta-\frac{q^2}{2m_\chi}, \label{Ee95def} \end{align}\tag{27}\] where \(\theta\) is the angle between \(\boldsymbol{q}\) and \(\boldsymbol{v}\). Solving for the incident speed gives \[\begin{align} v =\frac{1}{\cos\theta} \left( \frac{E_e}{q}+\frac{q}{2m_\chi} \right) = {v_\text{min}\over\cos\theta}. \end{align}\] For fixed \(q\) and \(E_e\), the smallest kinematically allowed speed is obtained for \(\cos\theta=1\).

For fixed \(v_\text{min}\) and \(E_e\), the equation \[\begin{align} E_e=qv_\text{min}-\frac{q^2}{2m_\chi} \end{align}\] has the two solutions \[\begin{align} q_-(v_\text{min},E_e) &=m_\chi v_\text{min} -\sqrt{m_\chi^2v_\text{min}^2-2m_\chi E_e}, \nonumber\\\; q_+(v_\text{min},E_e) &=m_\chi v_\text{min} +\sqrt{m_\chi^2v_\text{min}^2-2m_\chi E_e}. \label{q43q-} \end{align}\tag{28}\] The function \(v_\text{min}(q, E_e)\) has a unique minimum \(v_*=v_\text{min}(q_*,E_e) =\sqrt{{2E_e}/{m_\chi}}\), at \(q_* = \sqrt{2m_\chi E_e}\). Real solutions exist when \(v_{\rm min} \geq v_{\ast}\).

Consequently, for any \(v_\text{min}\) above the minimum value there are two momenta \(q\), namely \(q_-\) and \(q_+\), satisfying \(q_-\!\left(v_\text{min}\right) < q_* < q_+\!\left(v_\text{min}\right)\). Thus, to change the integration variable from \(q\) to \(v_\text{min}\), an integral over a range \(\left[q_i,q_f\right]\) containing \(q_\ast\) must be divided into two branch intervals on which the mapping is one-to-one. Their derivatives are \[\begin{align} \frac{\partial q_-}{\partial v_\text{min}} &= m_\chi- \frac{m_\chi^2v_\text{min}}{\sqrt{m_\chi^2v_\text{min}^2-2m_\chi E_e}} <0, \nonumber\\\; \frac{\partial q_+}{\partial v_\text{min}} &= m_\chi+ \frac{m_\chi^2v_\text{min}}{\sqrt{m_\chi^2v_\text{min}^2-2m_\chi E_e}} > 0. \end{align}\]

Defining \(v_i = v_\text{min}\!\left(q_i\right)\) and \(v_f = v_\text{min}\!\left(q_f\right)\) and reversing the limits on the lower momentum branch gives \[\begin{align} I &= \int_{q_i}^{q_f}\text{d}q F\!\left(q\right) = \int_{q_i}^{q_*}\text{d}q F\!\left(q\right) + \int_{q_*}^{q_f}\text{d}q F\!\left(q\right)\nonumber\\ &= \int_{v_*}^{v_i}\text{d}v_\text{min}{\left|\partial q_-\over\partial v_\text{min}\right| F\!\left(q_-\!\left(v_\text{min}\right)\right)}\nonumber\\ &\quad+ \int_{v_*}^{v_f}\text{d}v_\text{min}{\left|\partial q_+\over\partial v_\text{min}\right| F\!\left(q_+\!\left(v_\text{min}\right)\right)}~. \end{align}\] Thus, the absolute Jacobians account for the opposite orientations of the two branches.

For the scattering rate integrals considered in this work the original momentum range is \(q\in\left[0,\infty\right]\), for fixed \(E_e>0\) both transformed upper limits are formally infinite. Because the halo function \(\widetilde{\eta}(v_\text{min})\) vanishes above the maximum speed \(v_{\text{max}}\) supported by the detector frame DM distribution both integrals can instead be truncated at \(v_{\text{max}}\). For the SHM with fixed Earth speed \(v_e\), \(v_i =v_f=v_{\text{max}}=v_{\text{esc}}+v_e\) (which defines \(q_i = q_-(v_{\text{max}})\) and \(q_f = q_+(v_{\text{max}})\)).

9 Derivation of Event Rate↩︎

This derivation follows closely Refs. [69], [70]. Starting from the DM interaction rate defined in Eqs. 1 4 , imposing energy conservation for the deposited energy \(E\), we obtain

Using the two momentum transfer branches derived in App. 8, the \(q\) integral can then be transformed into the \(v_\text{min}\) integral given in Eq. 14 .

10 Truncated Maxwellian Halo Function↩︎

Each benchmark model of Sec. 6.1 makes use of the truncated Maxwellian velocity distribution of the form of Eq. 17 , differing only in the parameters \((v_0, v_e, v_\text{esc})\) and the density prefactor. Here, we state the corresponding closed form halo function expression, obtained by substituting Eq. 17 into the definition Eq. 9 and carrying out the angular integration.

For parameters \((v_0, v_e, v_\text{esc})\) and density \(\rho\),

where the truncated-Maxwellian normalization factor is \[\begin{align} K(v_0, v_\text{esc}) = \pi^{3/2} v_0^3 \left[ \text{erf}\!\left(\frac{v_\text{esc}}{v_0}\right) -\frac{2}{\sqrt{\pi}}~\frac{v_\text{esc}}{v_0}~ e^{-v_\text{esc}^2/v_0^2} \right]. \label{Kfactor} \end{align}\tag{30}\] The benchmark halo model functions then follow directly. For the Earth-bound benchmark for which \(v_e^{\rm EB}=0\) Eq. 29 is understood in the smooth \(v_e\to0\) limit.

References↩︎