Signatures and Detection Prospects for sub-GeV Dark Matter with Superfluid Helium

We explore the possibility of using superfluid helium for direct detection of sub-GeV dark matter (DM). We discuss the relevant phenomenology resulting from the scattering of an incident dark matter particle on a Helium nucleus. Rather than directly exciting quasi-particles, DM in this mass range will interact with a single He atom, triggering an atomic cascade which eventually also includes emission and thermalization of quasi-particles. We present in detail the analytical framework needed for modeling these processes and determining the resulting flux of quasi-particles. We propose a novel method for detecting this flux with modern force-sensitive devices, such as nanoelectro-mechanical system (NEMS) oscillators, and derive the sensitivity projections for a generic sub-GeV DM detection experiment using such sensors.


Introduction
Although the existence of dark matter (DM) has been definitively supported by gravitational evidence [1][2][3], its nature is still one of the biggest mysteries in modern physics.Recent theoretical developments [4,5] and tighter exclusion limits from both direct detection experiments [6][7][8][9][10][11][12] and the Large Hadron Collider (LHC) at CERN [13] are increasingly pointing towards DM being lighter than first thought, with mass less than ∼1 GeV.This motivates the design of new direct detection experiments which specifically target the sub-GeV and sub-MeV DM mass range [14][15][16].
Direct detection experiments achieve the greatest sensitivity when the fundamental excitation of the target material corresponds to the recoil energy or momentum scale from the scattering dark matter.Therefore, superfluid 4 He, with a spectrum of quasi-particle excitations with momenta keV [17,18], becomes an excellent target material candidate for detecting very light dark matter [19][20][21][22][23][24][25][26][27].Although superfluid 4 He has been considered as a target for neutrino and dark matter studies since the 1980s [28][29][30][31][32][33][34][35][36][37], its experimental application developed swiftly in recent years as the particle physics community refocused its interest on the search for light dark matter candidates [5,14,15,[38][39][40][41][42][43][44].One effort in applying superfluid 4 He is the search for light Weakly Interacting Massive Particles (WIMPs) of mass below 10 GeV [45][46][47].Due to the kinematics, the minimum velocity of light WIMPs required for elastic nuclear recoil is fairly low for Helium compared to heavier materials like Xenon and Germanium, and this fact could offer some additional sensitivity beyond that of the Xenon experiment.Other efforts have focused on the search for light dark matter candidates with masses below MeV [20][21][22][23][24].They utilize the fact that the fundamental excitation (quasi-particle) spectrum of the superfluid 4 He matches the recoil energy scale or momentum scale of the scattering dark matter.Therefore, each dark matter scattering event leads to the production of one or two phonon quasi-particles in the superfluid, for which the event rates and cross sections can be derived analytically.However, the practical detection of single excitations in the superfluid remains an experimental challenge.
In lieu of these previous proposals, in a previous theoretical work [48] we studied the possibility to use superfluid 4 He to detect DM in the sub-GeV mass range, i.e., for a range of DM masses between 1 MeV and 1 GeV. 1 The DM of our Milky way local group has a typical velocity of magnitude 10 −3 c [50 -54].Therefore, a DM particle with mass above 1 MeV has a de Broglie wavelength smaller than ∼ O(1) Å, the distance between helium atoms in the superfluid.As a result, the DM initially scatters with a single helium atom rather than producing multiple phonon quasi-particles.As we shall elaborate, this mass range conveniently produces a predominantly neutral atomic cascade, which has the advantages of (a) increasing the yield of quasi-particles and improving the sensitivity compared to that in the sub-MeV range; and (b) simplifying the modeling and projections of the experimental reach.
The basic physics processes are illustrated in the pentaptych of figure 1.In the first panel, the recoiling helium atom inherits the ∼ MeV momentum scale of the DM.It will proceed and scatter against other helium atoms, which themselves will scatter in turn, and so on, producing an avalanche of atoms (second panel).At the end of this atomic cascade, the initial recoil momentum has been distributed over a large number of daughter atoms moving isotropically with respect to the initial scattering location.These atoms have a de Broglie wavelength comparable or larger than the inter-atomic distance of ∼ O(1) Å.At this level, the 2-2 elastic scattering becomes subdominant and each low momentum atom interacts with the superfluid background and radiates quasi-particle excitations (including but not limited to phonons), as shown in the third panel.In our previous theoretical work, we focused on constructing an effective field theory (EFT) to describe this process [48].Here we build a numerical simulation which explicitly models the (relevant processes in the) atomic cascade.We will find that the quasi-particles are produced in a small region ∼ O(1) nm around the DM impact, and subsequently thermalize (fourth panel).Using the thermal distribution of quasi-particles we can trace the transport of momentum through the superfluid and derive the momentum imparted to a generic nanoelectro-mechanical system (NEMS) sensor suspended within (fifth panel).This allows us to derive the sensitivity projections for a generic sub-GeV DM detection experiment using a NEMS sensor.
Figure 1: Schematic illustration of a dark matter particle with mass between 1 MeV and 1 GeV scattering off superfluid 4 He.The de Broglie wavelength λ DM of the incoming DM particle is smaller than the inter-atomic spacing a, so that it scatters with a helium atom at rest (green dot).The dashed green lines correspond to scattered helium atoms freely propagating until being scattered by other helium atoms (green dots).This cascade process becomes subdominant when the de Broglie wavelengths λ He of the scattered helium atoms are comparable or larger than the atomic spacing, i.e., when λ He > a. Quasi-particles emitted by the He atoms are plotted with red lines.The quasi-particles rapidly interact with each other via various channels (depicted by the Feynman diagrams).Finally, a flux of quasi-particles with a Bose-Einstein distribution reaches the detector (yellow square).
The organization of our paper is as follows.In section 2, we review the elastic scattering between a DM particle and an initial helium atom.The calculation produces the helium recoil rate profiles for DM of different masses.In section 2.2, we review other scattering processes that may happen in the superfluid, and explain that they are subdominant for the sub-GeV regime.In section 3, we review the well known quantum mechanics treatment of helium-helium atomic cascade.In section 3.2, we review our previous theory proposal's result on the Lagrangian construction of helium atoms emitting quasi-particles.The calculation provides a quasi-particle production rate for a helium atom of an arbitrary momentum.In section 4.1, we review the parameters of quasi-particles and show that they are thermalized by the time all quasi-particles are radiated from the slow helium atoms.In section 4.2, we determine the temperature of the thermalized quasi-particle system by two different methods -an analytical approximation or a Monte Carlo simulation.In section 5, we calculate the momentum signal on the oscillator sensor in a realistic spatial configuration.Finally, we present the experimental constraints on the DM-Nucleus coupling strength.The five sections are chronologically ordered, with each section feeding a particle profile to the next section (see fig. 1).In section 6, we conclude our findings.
2 Dark matter scattering off Helium atoms

DM elastic scattering
We begin by considering the DM scattering rate off helium atoms via elastic nuclear recoil, which is the dominant process in the sub-GeV DM mass range.As shown in fig. 1, the superfluid Helium response to a dark matter scattering event depends only on the momentum of the initial recoiling He atom and is agnostic to the microscopic physics of the dark matter sector.For the purposes of our analysis, in this section we shall introduce a simple toy model of dark matter interactions which will be used to derive experimental sensitivity limits in section 5. Suppose that a fermionic dark matter particle χ couples to the nucleon n (we assume equal couplings to neutrons and protons) via a heavy scalar mediator particle φ, The differential DM-Helium cross section is where A = 4 is the atomic number of Helium.To arrive at the last result, we have taken the massive mediator limit m φ q max and approximated the nuclear form factor as F (q 2 ) → 1, which is justified by the fact that the inverse momentum of a sub-GeV mass DM, /q, is much larger than the nucleus size ∼ fm.
The differential nuclear recoil rate dΓ/dq 2 is linearly proportional to the number of target helium atoms N He , the local DM density ρ = 0.3 GeV/cm 3 , and depends on the DM velocity distribution f χ (v) and the differential cross section (eq.(2.2)), where in the second line we have introduced the DM-nucleon reduced mass, µ χn , and the DM-nucleon cross section σ χn = . The DM velocity distribution f χ (v) in our local frame is taken as the Maxwell-Boltzmann distribution with an escape velocity cutoff v esc , where the normalization factor The Heaviside function in eq.(2.4) constraints the DM velocity in the galactic frame to be below the escape velocity.Here we choose the velocity dispersion  The recoil momentum of the Helium atom is the essential quantity which governs the spectrum of the produced quasi-particles.By using eqns.(2.2-2.4),we obtain the differential Helium recoil rate shown in fig. 2. We note that heavier dark matter particles imply harder recoil spectra.

Other DM scattering signatures: direct production of quasi-particles and inelastic scattering
In addition to elastic scattering off individual helium atoms, there are several other DMinduced processes, which are subdominant in the MeV-GeV DM mass range targeted in this paper, but may become dominant for DM masses either below MeV or above GeV.In the sub-MeV mass range, the DM de Broglie wavelength is larger than several Å, spanning several helium atoms in the liquid.Such a light DM particle could directly produce quasiparticle excitations via coherent scatterings.At higher masses above a GeV, the exchange momentum is large enough to trigger excitation and ionization of helium.
Coherent quasi-particle production.The general scattering rate involving both coherent and incoherent quasi-particle production can be derived using the many-body quantization method [20,48]   Here we take a dark matter-nucleon cross section of 10 −40 cm 2 and 1 kg of target superfluid.The solid green line represents the event rate of single and multiple quasi-particle production with a net energy threshold of 1.2 meV.The nuclear recoil assumes no nuclear recoil below 4.5 keV.The dot-dashed orange data is taken from [22], where the 2-phonon detection threshold is at an energy of 1 meV.For a discussion of these threshold values, see refs.[20,22].
where the form factor S(q, ω) is the Dynamic Structure Function (DSF) of the superfluid.At low momentum, q ∼ keV, the DSF is measured from the neutron energy loss in neutron scattering experiments [55].At high momentum, q keV, the scattering is mainly incoherent, which is initiated by nuclear recoil and produces quasi-particles by the subsequent helium radiation.The integration of the incoherent part of the DSF dωS inc (q, ω) = 1 [56], so that the quasi-particle production rate is identical to the nuclear recoil rate in eq.(2.3), as expected.In principle, the DSF S(q, ω) can be applied to study DM scattering at any momentum, but at high momentum S(q, ω) is unknown and it is practical to use the nuclear recoil formula eq. ( 2.3).The coherent scattering becomes dominant at low momentum, where the measured or simulated DSF [55,57] could be employed to evaluate the relevant quasi-particle production rate.The production rate of phonon quasiparticles can also be derived using either the impurity method [48,58] or effective field theory [22,59,60].
To compare elastic scattering with coherent emission of quasi-particles, we show their respective event rates in fig. 3. The nuclear recoil rate (blue dashed line) is dominant for masses above ∼ 1 MeV, while coherent emission (green solid line) dominates for sub-MeV mass DM.The nuclear scattering rate is derived from eq. (2.3) by integrating over the momentum q 2 starting from q c = 4.5 keV, assuming that below q c DM will scatter Figure 4: A summary plot of helium atom elastic scattering, helium ion charge exchange scattering, and helium ion ionization cross sections [45], which interpolate the results in .The curves represent helium elastic scattering (blue solid), He 1+ to He 0 (He 2+ to He 0 , He 2+ to He 1+ , He 0 to He 1+ ) charge exchange scattering (dashed lines), and He 0 ionization (brown dot-dashed).The cross sections of other possible channels are lower than 10 −3 Å2 .
coherently with the superfluid.The coherent scattering rate is derived from eq. (2.5) by integrating over q 2 up to q c = 4.5 keV, with the DSF data given in [57].The DSF S(q, ω) peaks at ω ∼ meV and decreases quickly at high energy, because producing a large number of quasi-particles from a single coherent scattering is suppressed by couplings and phase space.The exchanged momentum in coherent scattering is restricted up to q c , while the momentum of nuclear recoil is q m χ v. Then for masses larger than MeV the coherent rate scales as 1/m 3 χ , while the incoherent rate goes only like 1/m χ and is therefore dominant at high mass.Furthermore, in each event, the coherent scattering generates O(1) quasi-particles, but the incoherent scattering could produce many more quasi-particles, depending on the DM mass.
The recoil energy is dissipated in the following channels: (1) ionized atom conversion into neutral atoms, (2) neutral atom cascade, (3) quasi-particle production by atoms, (4) decay of excited atoms into IR photon and singlet and triplet dimer excimers A 1 Σ + u and a 3 Σ + u .Existing efforts in the literature [19,45,46,88,89] have led to a clear estimation of the partition of recoil energy among these interaction channels.Essentially, quasi-particle production is the only relevant process at recoil energies below 20 eV, and it accounts for the dominant fraction all the way up to 100 keV.

Helium atom cascade
In this section we will discuss the theoretical description and the details of our simulation of the neutral atomic cascade.We simplify the DM-helium scattering model by considering only the neutral atoms portion of the cascade.This is justified by the numerical analysis in [45] and the summary in section 2 showing that the primary recoiling helium atom/ion will quickly produce a neutral atomic cascade.

Elastic scattering of Helium atoms
A recoiling helium atom with momentum in the keV to MeV range has de Broglie wavelength smaller than the inter-atomic distance in the superfluid.Therefore, the initial scattered atom triggers a cascade, i.e., a series of 2 → 2 scatterings among helium atoms within the fluid.Because the kinetic energy is comparable to the helium atomic potential [90][91][92], the 2 → 2 atomic scattering at this energy scale is non-perturbative.Nonetheless, given a potential function, we may numerically solve the Schrodinger equation using Partial Wave Expansion.Consider the wavefunction Ψ(r) as a function of the displacement vector r between two helium atoms.The initial and final asymptotic boundary conditions constrain the wavefunction as follows: The first term is an initial plane wave, the second term is an outgoing spherical wave weighted with a scattering amplitude f Ψ (θ, k), and k is the reduced momentum in the center-of-mass frame (equal to one half of the incoming atom momentum in the lab frame).Expanding eq.(3.1) in terms of Legendre Polynomials P l , the Schrodinger equation of the system can be solved numerically, and the scattering amplitude f Ψ (θ, k) depends on the phase shift δ l (k), where l is the angular index of the Legendre expansion.In [48] we calculated the scattering cross section using phase variation [93] and WKB approximation [94,95].We use the total scattering rate for the estimation in the next section of the helium cascade time which is related to the final shower size and quasi-particle number density.The differential cross section is incorporated in the simulation of the atomic cascade, and has the form For the convenience of comparison with the next section, in fig. 5 we plot the experimental total cross section using data from [91] (red dashed curve).

Quasi-particle emission from Helium atoms
The elastic scatterings keep dissipating the energy to more and more helium atoms, and decreasing their momenta.When a helium atom's momentum drops below 10 keV, the atom's de Broglie wavelength becomes comparable to the inter-atomic distance of the superfluid.The moving atom then collectively scatters against the surrounding atoms in the fluid and produces quasi-particles.
Unlike the case of sub-MeV dark matter directly producing quasi-particles [20][21][22]24], the helium emission process is non-perturbative.In [48] we proposed an effective U (1) current-current coupling between a helium atom and the superfluid: Λ is the cutoff scale of the superfluid EFT.We can estimate Λ using the inverse of the atomic space, Λ ∼ O(1) keV.For the phonon, the cutoff is related to the energy density of the superfluid ρ and the sound speed of phonon, Λ = (ρ c s ) 1/4 0.83 keV.J 0 is the number density operator of the superfluid, and its normalized matrix element is the Dynamic Structure Function (DSF) [20,23,96].The general form of the total emission rate of multiple quasi-particles is as follows: where ρ is the superfluid mass density, v He is the initial helium velocity, and S(k, ω) is the DSF.
For the purpose of comparing the elastic atomic scattering cross-sections with that for quasi-particle emission, in fig. 5 we show the emission rate for single quasi-particles.This is because our numerical integration of the phase space for multiple production ∼ d 3 k S(k, ω) and of the phase space for single production ∼ d 3 k S(k) δ(ω − ω(k)) with current available data shows that single quasi-particle production is dominant.In fig. 5 we show results for several combinations of the coupling parameters: λ 1 = 4π, λ 2 = 0 (blue), λ 1 = 0, λ 2 = 4π (orange), and λ 1 = 4π, λ 2 = 4π (green).Although the values of λ 1 and λ 2 are unknown, in this region the system is strongly coupled, so we take the value of 4π.In the simulation, we choose λ 1 = 4π, λ 2 = 0, but when we consider the energy/momentum conservation and thermalization after the quasi-particle production, the final quasi-particle flux will not be sensitive to the specific couplings due to the thermalization.The physical ) Elastic Exp.data

Figure 5:
The rate of helium emitting one quasi-particle (green, orange, and blue lines) and the experimental rate of elastic atomic scattering (red dashed line).The helium atom momentum is the lab frame momentum, consistent in all curves.The elastic atomic scattering rate is converted from the experimental scattering cross section given in [91].We test three combinations for the coupling parameters in eq.(3.5), namely, λ 1 = 4π, λ 2 = 0 (blue), λ 1 = 0, λ 2 = 4π (orange), and λ 1 = 4π, λ 2 = 4π (green).The orange and green lines match the expectation that at low momentum helium predominantly radiates quasi-particles.
expectation of low energy scattering originates from the two requirements: (1) helium atoms predominantly radiate quasi-particles rather than elastically scatter with each other in the superfluid; (2) quasi-particle modes only exist below ∼ 7 keV.The first requirement shows that the vector coupling in eq.(3.4) must exist.The second requirement restricts the momentum integration to an upper cutoff ∼ 7 keV.Below the cutoff, the emission rate increases with the momentum of the atom because of the phase space enhancement.Beyond the cutoff, the emission rate decreases because the phase space integration in eq.(3.5) has reached its maximum, leaving powers of helium atom momentum only in the denominator of the result.The curves in fig. 5 thus have a cusp at 7 keV because of this different behavior below and above the cutoff.

Monte Carlo simulation of the Helium cascade and radiation of quasiparticles
In this section, we will present our simulation results about the developing helium cascade and the radiation of quasi-particles from fast moving helium atoms.We do not include the subsequent quasi-particle decays and quasi-particle self-interactions, whose treatment is postponed for the next section.
Because of the large number of daughter particles which must necessarily be produced to conserve both energy and momentum, little information about the overall final structure of the shower may be gleaned from a direct analytical treatment.However, because the individual interactions of the shower constituents are quantum mechanical and probabilistic, the problem is amenable to a Monte Carlo approach.In particular the distance that each helium atom travels between interactions, the type of interaction it experiences, and the subsequent evolution of its daughter particles are all properly described by probability distributions (all necessary results were derived and/or collected in [48]), which we exploit to generate ensembles of simulated events.
The momentum-dependent cross section σ el of elastic atomic helium scattering is known from experiment [91] (see red line in fig.5).On the other hand, the rate of inelastic emission of quasi-particles has been computed in [48] and is given by eq.(3.5).This rate can be cast as an inelastic "cross section" according to the heuristic where n He = ρ/m He denotes the number density of helium atoms in the superfluid.The total cross section of these processes defines a mean free path and thus a probability distribution that describes the distance an energetic helium atom is expected to travel before interacting with the superfluid to produce either another energetic helium atom or a quasi-particle.The respective probabilities for each type of daughter particle are given by Note that the cross sections and consequently the length scale of the probability distribution are all functions of the helium momentum.At this point we have all the ingredients needed to describe the structure of the Monte Carlo simulation: each simulated event begins as a single helium atom at the origin with its momentum oriented along the z-axis.
Using this momentum we evaluate the cross sections of both processes and sample the resulting distribution to determine what type of interaction the helium atom experiences and where this interaction occurs.The momenta of daughter particles are sampled from the appropriate differential rates provided in sections 3.1 and 3.2 above, and the cross sections of those daughter particles are evaluated anew.In this way the simulation proceeds recursively generating new daughter particles and tracking their trajectories between collisions.This process is depicted schematically in the second and third panels of fig. 1.
In fig.6 we show typical quasi-particle spectra from our simulation.In the first (second, third, fourth) row of panels we show distributions of quasi-particle momenta (energies, velocities, cos θ), for an initial He atom momentum of 50 keV (left column) and 500 keV (right column).The top panels show that the number of phonons is significantly lower  (compared to the number of rotons) due to the phase space suppression.Nonetheless, quasiparticle decays to softer modes, not included in our simulation, will eventually populate that region.The panels in the second row reveal that the energy spectrum starts with a peak at the gap energy ∆ = 0.75 meV which corresponds to the local minimum of the roton dispersion.This is also reflected in the velocity graphs on the third row, where the peak around zero velocity is composed of slow rotons and maxons, as well as slow quasiparticle modes above ∼ 5 keV.The plots in the last row demonstrate that the resulting distribution is almost isotropic.The discussion in section 4.1 below will show that all these quasi-particles will become thermalized.Therefore, instead of using the generated quasiparticles from fig. 6, in the later sections we shall sample the quasi-particle spectrum from a Bose-Einstein distribution.

Effects of interactions among quasi-particles
In this section we investigate the effects of interactions between the quasi-particles produced at the end of the cascade.These interactions may thermalize the ensemble of quasiparticles.In order to determine whether this actually occurs, we first do a back-of-theiphone estimate comparing the quasi-particle interaction rate nσv and the quasi-particle production rate Γ inel , for a given recoil energy or momentum.After that, we present the theoretical solution of last scattering surface and the thermal temperature.More accurate results from our simulation are presented in section 4.3.

Theoretical overview
We consider the two well known types of excitations with analytic dispersion relationsphonons and rotons.The "phonon" refers to a quasi-particle with momentum below ∼ 1.2 keV; its dispersion is linear with a small cubic correction: where c s = 240 m/s is the sound speed, Λ = (ρc s ) 1/4 is the UV scale of the superfluid, and γ ∼ O( 1) is a dimensionless parameter such that γ/Λ2 = 0.27 Å2 .The "roton" is a quasi-particle of momentum ∼ p * = 3.84 keV at the local minimum of the dispersion curve.Its energy is parameterized as where m * 0.16m He is the effective roton mass.Several relevant interactions are well studied in the literature [18, 48, 58-60, 97, 98], including phonon decay, phonon 2-2 scattering, roton 2-2 scattering, and phonon-roton 2-2 scattering 2 .In Table 1 of [48], we listed the main results for these interaction cross sections.
With those ingredients, we are now ready to check for thermalization.If at the end of the atomic cascade, when all quasi-particles have just been produced, the quasi-particles have already experienced multiple interactions, we can safely claim that the quasi-particle system is thermalized.We perform a back-of-the-envelope estimation as follows.
From fig. 2, we know that the initial recoil momentum P ini of the helium atom triggering the cascade ranges from 1 to 10 3 keV.According to fig. 5, when the helium atom momentum drops to below ∼ 20 keV, the atom will predominantly start to radiate quasi-particles, thus not affecting the number of atoms in the cascade.Therefore, energy conservation implies that a recoil momentum of P ini keV will dissipate to P ini 20 keV 2 slower helium atoms.Each of those slow atoms in turn will radiate about 100 quasi-particles, assuming that the quasiparticles' energies are ∼ 1 meV.We take the typical scattering and radiation rate from fig. 5 as 10 13 s −1 .The radiation of all quasi-particles will take ∆t ∼ 100 × 10 −13 s = 10 −11 s, which is longer than the prior atomic cascade time because the number of atoms increase exponentially with time during a cascade.Therefore, we estimate the total time of cascade and radiation to be 10 −11 s.
During this time, the quasi-particles (with velocity O(100) m/s) may expand up to ∆R ∼ 100 m/s × 10 −11 s = 10 −9 m = 1 nm.We then estimate the the interaction rate nσv for quasi-particle self-interactions as We see that the corresponding timescale ranges from 10 −14 s (for P ini ∼ 1 MeV) to 10 −11 s (for P ini ∼ 20 keV).This timescale is smaller than the timescale of 10 −11 s for producing the quasi-particle shower.Therefore, quasi-particles dissipated from P ini 20 keV recoil momentum will interact with each other multiple times, i.e. become fully thermalized, by the time when all quasi-particles are produced.When considering recoil momenta 20 keV, only a few quasi-particles are produced at this scale, and they infrequently interact with each other.Thus we treat these quasi-particles as free-streaming from the beginning.

Quasi-particle thermalization
In the following, we present our procedure to derive the thermalized distribution of quasiparticles.First of all, we estimate the last scattering surface of the final quasi-particle "plasma" system, assuming an isotropic configuration distribution.Borrowing from the analogous concept in Cosmology [99,100], at the last scattering surface, the quasi-particle interaction rate equals the expansion rate.Modeling the isotropic quasi-particle "plasma" as a sphere, the last scattering radius R ls is determined by the radius for which the optical depth is unity: where N is the total number of quasi-particles, σ is their interaction rate, and V = 4πr 3 /3.The solution of the previous equation gives the last scattering radius R ls = 3N σ/8π.N is estimated by assuming that all final quasi-particle excitations have energy 1 meV: N = P 2 ini /2m He 1 meV .Plugging in the numerical values from the previous subsection, we find the magnitude of the last scattering radius R ls to be O(1) to O(10) nm, smaller than the distance between the sensors in realistic detectors.Therefore, the relevant distribution detected by the sensors is the thermalized spectrum of quasi-particles.The number of quasi-particles is not fixed, but can change due to a) decays of unstable phonons (with momenta below 1.2 keV) to other phonons; b) decays of other quasi-particles to stable quasi-particles with momenta between 1.2 and 4.6 keV [19,101]; and c) multiple quasiparticle interactions which can change the number of quasi-particle modes.Therefore, the literature treats the chemical potential of quasi-particles as zero.Then we can express the Bose-Einstein distribution of the quasi-particles as Using energy conservation, the initial recoil energy equals the total energy of the Bose-Einstein ensemble: where 4  3 πR 3 ls is the volume of the thermal system.The momentum integration runs from 0 to 4.6 keV because we assume only phonons and stable quasi-particles exist in the thermalized distribution.

Monte Carlo simulation and the thermal temperature
We can use our MC simulations described in section 3.3 to provide a numerical estimate for T from the previous subsection.The simulation provides averaged values of the different scattering cross sections between quasi-particles of all types.These cross-sections can be used in eq. ( 4.3) for the calculation of R ls , which can then be substituted in eq.(4.5) to solve for T .The result is shown as the blue solid line in fig.7, for which we have used the value σ = 8.64 keV −2 suggested by our simulations.As mentioned in section 4.1, the region of P ini below 20 keV is unnecessary for our simulation, since the recoiled helium atom will radiate quasi-particles without an atomic cascade, and the number of quasi-particles will be insufficient for thermalization.
One potential problem with our simulations is that the expression for the phonon-roton scattering cross section is only valid for phonons that have a much smaller momentum than the roton [48,59,60].Therefore, we cannot reliably sample the cross section between hard phonons and rotons.In order to get an idea of the effect from this theoretical uncertainty, we introduce a k-factor for the phonon-roton cross section and in fig.7 show results for k = 10 (red dashed line), k = 1 (green dashed line) and k = 0.1 (orange dashed line).The width of the shaded band enclosed between the orange and red dashed lines is indicative of the corresponding uncertainty on the derived thermal temperature T , which should be kept in mind when discussing projections for the experimental sensitivity below.The thermal temperature T predicted from eq. (4.5) as a function of the initial recoil momentum P ini .The blue solid line corresponds to taking σ = 8.64 keV −2 (the value of roton-roton average cross section provided by the simulation) when solving for the last scattering radius R ls appearing in eq.(4.5).The dashed lines in the range 50 to 200 keV are obtained by using the simulated average σ, with the phonon-roton scattering cross section multiplied by a k-factor of k = 1 (green), k = 0.1 (orange), or k = 10 (red).

Experimental signals
The previous analysis and simulations are applicable to various types of dark matter searches in superfluid helium.One such proposal is the HeRALD experiment [19], which tries to detect dark matter scattering events occurring in the bulk of the superfluid by sensing quantum evaporation of helium atoms from the superfluid surface.Here we study the detection prospects of mechanical oscillators like NEMS, which instead sense the force (or momentum deposit) due to quasi-particles.For a given momentum threshold and detector size, we deduce the potential to discover the dark matter signal using mechanical oscillators and compare with the HeRALD experiment, leaving the detailed experimental design and study of the backgrounds for a future work [102].
Detector setup.In order to study the quasi-particle flux onto a generic device, we consider a simple detector model as follows.We assume a large supply of 2D squareshaped sensors of area a × a, which are placed on a 3D grid inside a superfluid target, with an inter-sensor distance of d.The distance between the event and the closest target sensor is up to half of a diagonal √ 3d/2.Each sensor at the cell vertex first receives the signal from its own octic cube.Therefore, we simplify the geometry and study the events happening within a cube of volume d 3 with a sensor at the center (see fig. 8).We assume that a single sensor is a square of size a × a, and the sensors are arranged in a 3D lattice with spacing d.The event can happen at any location in the lattice, but it will be first detected by the closest sensor.For one given sensor, we thus calculate the signal from events within a cube of volume d 3 .See section 5 for details.
We also model the detection as an event free streaming to the sensor of area a 2 from a distance r.Assuming the equilibrium quasi-particle ensemble to be isotropic (which is confirmed in our simulations), the probability of a single given particle reaching the surface is given by the ratio between its solid angle to the detector Ω(θ, φ, r, a) and the total solid angle 4π.The solid angle Ω(θ, φ, r, a) belongs to a leaning pyramid configuration and is solved numerically using Mathematica build-in commands.If P ini > 20 keV, the quasi-particles reach thermalization.Then we calculate the differential profile, dN/dp, the number of quasi-particles with given momentum p, according to a thermal distribution: where V is the last scattering volume estimated from eq. ( 4.3), while T is the thermal temperature discussed in fig. 7. Results for a few representative values of P ini are shown in fig.9.For P ini < 20 keV, the quasi-particles are free-streaming from their point of origin.
We estimate the number of events using energy conservation.According to the simulation results in section 3.3, for the estimation of the quasi-particle flux we further assume that most of the quasi-particles are rotons.The spectrum of differential particle number dN/dp as a function of quasiparticle momentum p.The thermalized quasi-particle systems are generated by a recoil momentum of 20 keV (blue), 100 keV (orange), and 1000 keV (green) respectively.The total number of quasi-particles is the area under the curves.We assume only phonons and stable quasi-particles exist in the thermalized distribution.
O(1) keV threshold.Sensors of momentum threshold ∼ O(1) keV are capable of detecting single quasi-particles.Taking 1 keV threshold as an example, as long as at least one quasi-particle reaches the sensor, the event can be counted as detectable.There is a higher probability of detection when the particles are streaming from a closer distance and a normal angle from the sensor.We thus extract this probability between the event and the nearest sensor: Here p min is the threshold momentum of the sensor, and N (p min ) is the total number of quasi-particles with momentum above threshold, which is integrated from eq. (5.1).We plot an example density plot for this probability function in the left panel of fig.10.The number of quasi-particles arriving at the sensor decreases when the events happen further from the sensor and away from its normal direction.Since the DM-Helium scattering rate eq.( 2.3) is independent to the coordinate space, we average the probability P (θ, φ, a, r) over the whole sensor cell and call it the Spatial Efficiency (SE): where radial distance r = x 2 + y 2 + z 2 , the zenith angle θ = tan −1 ( x 2 + y 2 /z), and the azimuth angle φ = tan −1 (y/x).The spatial efficiency eventually depends on the recoil  The initial recoil momentum from the DM is 100 keV, the sensor area is 100 × 100 µm 2 .
Here we vary the zenith angle θ along the horizontal axis, while setting the azimuth angle φ = 0.The probability increases as the event radial distance decreases and as the direction gets closer to the normal (θ = 0).momentum P ini , the sensor distance d, and the sensor area a 2 .We show the spatial efficiency function eq. ( 5.3) results for O(1) keV threshold scenario in fig.11a.Comparing the three curves, we notice that the spatial efficiency is generally proportional to the area a 2 of the sensor.This is because the probability (5.2) takes the limit N (p min )Ω(θ, φ, a, r)/4π at a large particle number N (p min ) 1.In addition to the nearest sensor, the other sensors further away will also receive some signal, which we roughly estimate to be on the order of 40% (for a threshold of 1 keV), but to be conservative, this correction will not be included in our sensitivity projections below.O(10) keV threshold and beyond.Sensors of momentum threshold ∼ O (10) keV are incapable of detecting a single quasi-particle.Taking 10 keV threshold as an example, we estimate that it takes ≥ 3 quasi-particles reaching the sensor to trigger a detectable event.The probability (5.2) can thus be generalized into a cumulative Binomial distribution: Here the number N is the total number of quasi-particles without the cap of p min .Similar formalism can be applied to O(100) keV threshold.For example, we estimate that at least 25 quasi-particles are required to trigger a 100 keV threshold sensor.Therefore, the  cumulative sum in eq. ( 5.4) is up to m = 24.The total number of quasi-particles N must exceed the minimum quasi-particle number that triggers the threshold, which is the maximum value of m in this equation.For example, there are 3 quasi-particles produced below 5 keV recoil momentum (no matter whether thermalized or not), thus the probability and the spatial efficiency vanishes.
Despite these changes, the spatial efficiency function follows the same formalism (5.3).As shown in fig.11b, the spatial efficiency in this scenario is about three orders of magnitude smaller than the previous case.Moreover, the spatial efficiency no longer scales with the area but with a 3 .It is therefore unnecessary to compute the contribution from the sensors further away because they are unable to compensate for the loss of sensitivity.
Sensitivity projection.DM with a given mass could generate a recoil momentum of any value with the rate distribution shown in fig. 2. The momentum threshold and other parameters of the detector determine the probability of an event at location (θ, φ, r) being detectable.Mathematically, the integration of dΓ/dq 2 from fig. 2 is thus weighted by the curves in the fig.11 when we calculate the total event rate.Although the spatial efficiency is a continuous function, the total detectable event rate vanishes when the total number of quasi-particles is less than the minimum requirement to trigger the threshold, as mentioned above.
In fig.12, we plot the estimated reach in DM-Nucleus recoil cross section σ χn for several detector configurations, based on 90 percent confidence level (2.3 events) per exposure time per target mass [21,[103][104][105][106][107], assuming zero background.For a one year exposure with one kilogram of target, the constraint on the DM-Helium scattering rate is: Figure 12: The solid lines show 95% C.L. exclusion regions for our proposed experiment for different detector parameters (assuming zero background).The legends label the three detector parameters in the following order: the exposure time and target mass, the area of the sensor, and the momentum threshold.The distance between two sensors is set as 1 cm.We also include the HeRALD results [19] (dashed lines) which consider the solar neutrino background.The parameters for the HeRALD results are: 1 kg-day exposure with 40 eV energy threshold (blue), 1 kg-yr exposure with 10 eV energy threshold (orange), and 10 kg-yr exposure with 0.1 eV energy threshold (green).
where q min is the cutoff by N > 1 in eq. ( 5.2) or N > m max in eq.(5.4), and q max is the recoil momentum upper limit from classical kinematics.The variation of exposure includes per day per kilogram, per year per kilogram, and per year per 10 kilogram, with the expected number of events scaling accordingly.A combination of smaller distance between sensors, larger area of the sensor, and longer exposure time allows to probe smaller cross sections.

Conclusions
We have proposed a new method of detecting sub-GeV DM using the spectrum of the quasi-particle excitations in superfluid 4 He generated as a result of a DM collision with a He nucleus.The key idea is to leverage modern force-sensitive devices, such as NEMS os-cillators, to detect the momentum flux of the quasi-particles.With a kg-year exposure, we have demonstrated that this superfluid experiment can strongly constrain the DM-Nucleus interaction within the sub-GeV mass range.Our study is also the first to theoretically study (1) the production of quasi-particles from the resulting helium atom cascade; (2) the thermalization and decoupling of quasi-particles; and (3) the temperature of the thermalized quasi-particle system in a superfluid DM detector.The findings include: • The DM collision initiates a cascade of helium atoms which gradually lose momentum by elastic scattering, and later, once they reach O(10) keV, by radiation of quasiparticles.
• Quasi-particles generated by a recoil momentum P ini O(10) keV will thermalize as a result of their self-interactions.For P ini O(10) keV, the quasi-particles are free streaming from the beginning.
• The temperature of the thermalized quasi-particles system generated by sub-GeV DM is O(1) Kelvin.
These results lead to a prediction of the detectability of events in the vicinity of a NEMS sensor.With a 10 kg-yr exposure and a 1 keV momentum threshold, the ideal detector setup is able to push the DM-nucleon exclusion region down to about 10 −41 to 10 −43 cm 2 .

v 0 =
220 km/s, the velocity of the earth v e = 240 km/s, and the escape velocity v esc = 500 km/s.

Figure 3 :
Figure3: Event rates for the processes of nuclear elastic scattering (blue dashed line), coherent quasi-particle emission (green solid line), and double phonon emission (orange dot-dashed line).Here we take a dark matter-nucleon cross section of 10 −40 cm 2 and 1 kg of target superfluid.The solid green line represents the event rate of single and multiple quasi-particle production with a net energy threshold of 1.2 meV.The nuclear recoil assumes no nuclear recoil below 4.5 keV.The dot-dashed orange data is taken from[22], where the 2-phonon detection threshold is at an energy of 1 meV.For a discussion of these threshold values, see refs.[20,22].

Figure 6 :
Figure6: Typical quasi-particle spectra from our simulation.In the first (second, third, fourth) row of panels we show distributions of quasi-particle momenta (energies, velocities, cos θ), for an initial He atom momentum of 50 keV (left column) and 500 keV (right column).

Figure 7 :
Figure 7:The thermal temperature T predicted from eq. (4.5) as a function of the initial recoil momentum P ini .The blue solid line corresponds to taking σ = 8.64 keV −2 (the value of roton-roton average cross section provided by the simulation) when solving for the last scattering radius R ls appearing in eq.(4.5).The dashed lines in the range 50 to 200 keV are obtained by using the simulated average σ, with the phonon-roton scattering cross section multiplied by a k-factor of k = 1 (green), k = 0.1 (orange), or k = 10 (red).

Figure 8 :
Figure 8: Illustration of the geometry used to calculate the signal on the detector surface.We assume that a single sensor is a square of size a × a, and the sensors are arranged in a 3D lattice with spacing d.The event can happen at any location in the lattice, but it will be first detected by the closest sensor.For one given sensor, we thus calculate the signal from events within a cube of volume d 3 .See section 5 for details.
momentum (keV) Differential number dN/dp Generated by 20 keV recoil momentum Generated by 100 keV recoil momentum Generated by 1000 keV recoil momentum

Figure 9 :
Figure9: The spectrum of differential particle number dN/dp as a function of quasiparticle momentum p.The thermalized quasi-particle systems are generated by a recoil momentum of 20 keV (blue), 100 keV (orange), and 1000 keV (green) respectively.The total number of quasi-particles is the area under the curves.We assume only phonons and stable quasi-particles exist in the thermalized distribution.

Figure 10 :
Figure 10: Detection probability of two example sensors.(a) The probability of a 1 keV threshold sensor from eq. (5.2).The initial recoil momentum from the DM is 100 keV, the sensor area is 10 × 10 µm 2 .(b) The probability of a 10 keV threshold sensor from eq. (5.4).The initial recoil momentum from the DM is 100 keV, the sensor area is 100 × 100 µm 2 .Here we vary the zenith angle θ along the horizontal axis, while setting the azimuth angle φ = 0.The probability increases as the event radial distance decreases and as the direction gets closer to the normal (θ = 0).
Spatial efficiency of 10 keV threshold