{'fraction_non_alphanumeric': 0.037262758329818645, 'fraction_numerical': 0.01037536904259806, 'mean_word_length': 4.372266908349382, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Directional Detection of Light Dark Matter in Superconductors

Superconducting detectors have been proposed as outstanding targets for the direct detection of light dark matter scattering at masses as low as a keV. We study the prospects for directional detection of dark matter in isotropic superconducting targets from the angular distribution of excitations produced in the material. We find that dark matter scattering produces initial excitations with an anisotropic distribution, and further show that this directional information can be preserved as the initial excitations relax. Our results demonstrate that directional detection is possible for a wide range of dark matter masses, and pave the way for light dark matter discovery with bulk superconducting targets.

Superconducting detectors have been proposed as outstanding targets for the direct detection of light dark matter scattering at masses as low as a keV. We study the prospects for directional detection of dark matter in isotropic superconducting targets from the angular distribution of excitations produced in the material. We find that dark matter scattering produces initial excitations with an anisotropic distribution, and further show that this directional information can be preserved as the initial excitations relax. Our results demonstrate that directional detection is possible for a wide range of dark matter masses, and pave the way for light dark matter discovery with bulk superconducting targets.
The identity of dark matter (DM) remains one of the most pressing open questions in particle physics and cosmology. Contrary to decades of theoretical expectations, numerous experimental probes have found no conclusive evidence of DM at the weak scale, leading to renewed interest in models of DM at much lower scales. Many new ideas have recently been proposed to search for such light DM in the laboratory , and several of these novel direct detection exeperiments have already begun to probe significant parameter space [39][40][41][42][43][44][45][46][47][48][49][50][51][52]. Among the new ideas, superconducting targets stand out with the lowest possible thresholds, giving them sensitivity to the lowest DM masses through DM-electron interactions [5,6,24,53]. With superconducting energy gaps of O(meV), such detectors may eventually probe DM with mass as low as the keV scale, where cosmological constraints become significant [54][55][56].
Despite the impressive potential reach of superconducting targets, current projections assume that the detectors are insensitive to the direction of incoming DM. Directional detection has long been recognized as a powerful tool in DM experiments, including those in the keV-GeV regime [27,38]: due to the halo wind, the local DM distribution is not isotropic in the laboratory frame, leading to a characteristic modulation of the signal that can be used to reject backgrounds and confirm a discovery. If superconducting detectors can be made sensitive to the direction of the incoming DM, then such targets will offer exceptional promise for future experiments. Such a detector would be capable of making a definitive discovery of DM as light as a keV, should a signal be observed.
In this work, we show that even isotropic superconductors are capable of directional detection via the angular distribution of the excitations produced by DM scattering. For such a measurement to be viable, two key features are required. Firstly, the direction of the initial excitations produced by the DM interaction should be correlated with that of the incoming DM particle. Sec-ondly, the secondary excitations produced by the initial excitations as they down-convert in the material should exhibit directionality correlated with that of the initial excitations. As we will show, both features indeed occur in superconducting targets, paving the way for directional detection of keV-scale DM.
This Letter is organized as follows. We begin by considering the initial scattering of DM with electronic states of a superconductor into excited quasiparticles. Next we consider the down-conversion of these initial excitations in the material into secondary quasiparticles and phonons, treating the general case with a new numerical code. We then present our results for directionality, and end with a discussion of experimental prospects.
Throughout this work, we use the following notation: for a 3-vector q, we write q = |q|. We use angles with two subscripts to denote the relative angle between the two axes specified by those subscripts. All other angles are defined relative to the DM wind axis. The Fermi energy and momentum are denoted by E F and p F , respectively. We use natural units where = c = 1.

DARK MATTER SCATTERING
We begin by demonstrating the directionality of the initial excitations produced in a DM scattering event in a superconducting target. This requires reformulating the description of the scattering process in terms of the appropriate degrees of freedom in the 'BCS' vacuum of the superconductor [57]. The DM scattering rate in superconductors was originally computed by Ref. [6] considering only large energy deposits compared to the superconducting gap, and so the DM-detector interaction was described in terms of the DM χ scattering with individual electrons, |χ |e − −→ |χ |e − . While this description is appropriate for computing the rate, it is not suitable for studying the kinematics at small deposits. Here, the appropriate degrees of freedom are Bogoliubov quasiparticles (QPs) [58], which are electron-hole superpositions.
In this description, the DM excites the BCS vacuum by pair-producing QPs, as |χ |0 BCS −→ |χ |QP 1 , QP 2 . The total momentum of these two QPs is the momentum transfer q imparted by the DM scatter. The wave functions of the electrons in the BCS vacuum automatically account for Pauli blocking through a coherence factor, which has significant support only when one of the QPs is below the Fermi momentum and the other is above. In the Supplemental Material (SM), we show that for energy deposits much larger than the superconducting gap energy, the scattering rate becomes identical to that for scattering with individual electrons: the matrix element for QP pair production becomes that for electron scattering, the coherence factor becomes a Pauli blocking factor, and the two QPs become an electron-hole pair. We use the labels 1 and 2 to refer to the two initial QPs produced in the scattering process. All quantities are invariant under interchange of QPs 1 and 2, so we use the label i for statements that apply to either label.
The QPs have a dispersion relation different from that of free electrons, of the form where E p = p 2 /(2m * ) − E F is the Bloch energy relative to the Fermi surface and ∆ is half of the superconducting gap energy. Counterintuitively, the energy of a QP is minimized (E QP = ∆) when its momentum lies on the Fermi surface. The free-electron dispersion relation is recovered in the limit p QP p F , whereas the limit p QP p F gives the energy of a hole far below the Fermi surface. This non-trivial dispersion relation modifies the kinematics of DM scattering near the gap, and thus significantly influences directional correlations and downconversion. For energy deposits ω keV, the momenta p 1 and p 2 will be well inside the first Brillouin zone (BZ). For deposits ω keV, we expect Eq. (1) to receive band structure corrections near the edge of a BZ of order tens of eV or less, small compared to this keV scale, so Eq. (1) is a valid approximation for all the energy scales considered in this work.
The overall DM scattering rate is given by [53] where v χ = |v χ | is the magnitude of the DM velocity; q is the momentum transfer; ω q = q · v χ − q 2 /2m χ is the deposited energy; and BCS is the dielectric function of a superconductor in the BCS vacuum. In this work, we make the approximation Here L is the Lindhard form of the dielectric function for a normal metal [59], which accounts for the effects of in-medium screening and collec-tive modes in the normal metal phase [53]. We compute Im ( BCS ) in terms of the QP dispersion relation Eq. (1) and the BCS coherence factor, which accounts for neargap effects [60] and for Pauli blocking. Our approach interpolates between the approximate superconducting dielectric response near the gap and the normal-metal response far from the gap. A more complete treatment explicitly computing the dielectric function in the BCS vacuum lies beyond the scope of this work and will pursued elsewhere [61]. We take |V (q)| 2 = (g e g χ ) 2 q 2 + m 2 φ −2 , where g e and g χ are the couplings of the mediator to the electron and the DM, respectively. This is appropriate for any spin-independent interaction. Further details on DM interactions and the computation of Im ( BCS ) are given in the SM.
The rate in Eq. (2) can be differentiated to obtain the parameter distributions of the initial QPs prior to down-conversion. We begin with the joint distribution of cos θ χ1 , cos θ χ2 , E 1 , and E 2 , where θ χi is the angle between p χ and p i . While the distribution of the cos θ χi fully describes the directionality of the initial excitations, the energies of these excitations are also needed to determine the properties of secondary excitations. The four-parameter joint distribution f QP [p χ ](cos θ χ1 , cos θ χ2 , E 1 , E 2 ) at fixed p χ can be combined with a DM distribution f χ (p χ ) to obtain the full distribution f QP (cos θ 1 , cos θ 2 , E 1 , E 2 ), where θ i is the angle between the QP momentum and the DM wind axis. In our numerical results, we draw samples from this distribution to generate initial excitations.
In certain regimes, the angular distribution of excitations can be understood analytically by virtue of kinematical constraints. Conservation of energy yields a closed-form expression for cos θ qi , the cosine of the angle between QP i and the momentum transfer q, in terms of q and p i . Now consider small DM masses and light mediators, where small deposits are favored. In the limit of small deposits, to leading order in ω − 2∆, we have p i p F . Then conservation of momentum requires where v χ = |v χ |, which implies 0 cos θ qi 2∆/(v χ p F ). For aluminum (Al), with ∆ = 0.3 meV and E F = 11.7 eV, this leads to the condition 0 cos θ qi 10 −4 , so excitations produced near the gap are nearly orthogonal to the momentum transfer q. In turn, the direction of q is necessarily correlated with that of v χ , so the distribution of cos θ χi is peaked at zero.
On the other hand, consider DM interacting via a heavy mediator, for which large deposits are favored. In particular, for p χ ∼ p F , momentum transfers of order p F are possible, corresponding to cos θ qi ∼ 1. For example, for p χ = 2p F , if the DM is fully stopped and its energy shared equally between the two QPs, then cos θ qi is given uniquely by cos In Al, this solution corresponds to cos θ qi ≈ 0.93. Further, fully stopping the DM implies that cos θ χq = 1, so cos θ χi = cos θ qi . Thus, when p χ ∼ p F and large q is favored, the angular distribution can peak in the direction of the DM wind.
The marginal distribution of cos θ i is shown by the solid curves in the top panels of Fig. 1 for several DM masses in Al. We assume the Standard Halo Model (SHM), i.e., in the galactic frame, taking v 0 = 220 km/s, v esc = 550 km/s, and Earth velocity v E = 230 km/s relative to the galactic frame. For light DM or a light mediator, small energy transfers are favored, leading to a peak in the angular distribution orthogonal to the DM wind axis. For heavier DM and mediators, larger energy transfers lead to a forward-peaked distribution.

QUASIPARTICLE RELAXATION
The scattering process described above produces a single pair of QPs with a known total energy and net momentum, and we showed how the directions of these quasiparticles are related to the direction of the incoming DM. The second requirement for directional detection is that this directionality must be preserved after the initial excitations relax. Thus, we now study the down-conversion of the initial QP excitations through QP-phonon scattering.
Following Refs. [62][63][64], we model down-conversion as a repeating sequence of two distinct processes: first, energetic QPs relax by emission of phonons, and second, energetic phonons decay into QP pairs. Quasiparticle pair production eventually stops once all remaining phonons have energy below 2∆. We treat such phonons as ballistic, including them as part of the final state that is eventually read out by the detector. As shown in the SM, phonon emission is kinematically forbidden for very low-energy QPs, so the QPs eventually become ballistic as well. The down-conversion process is finished when all particles are ballistic. For other approaches to the relaxation of highly energetic QPs, see Refs. [65,66].
We study the impact of down-conversion on directionality by explicit simulation, implemented in a public code released together with this Letter. 1 We begin with an ensemble of initial excitations sampled from the distribution f QP , and then iterate the relaxation processes described above until all QPs and phonons are ballistic. Computing the momentum of the outgoing excitations after each relaxation process requires knowledge of the differential rate of these processes with respect to the kinematical parameters of the final state. We take the differential rates for phonon emission and QP pair production from Eqs. (16) and (27) of Ref. [67]. In each case, imposing conservation of energy and momentum using the dispersion relation of Eq. (1) gives the distribution of final-state angles. The differential rate of QP pair production by a phonon of energy ω ph is given by where E i is the energy of one QP and ω ph − E i is the energy of the other. This distribution is sharply peaked at E i ∼ ∆ and E i ∼ ω ph − ∆, corresponding to the case in which one of the two QPs receives most of the phonon's energy. Here, the dispersion relation of Eq. (1) implies that the QPs are produced nearly orthogonal to the axis of the initial phonon. Thus, for small deposits, the angular distributions of QPs and phonons in the final state will peak in orthogonal directions. We sample initial excitations and simulate the downconversion process for several DM masses in Al. The angular distributions of the resulting QPs are shown as the solid curves in the bottom panels of Fig. 1. The angular distribution of phonons exhibits weaker directionality and is discussed in the SM. We learn that for light mediators, the initial QPs have excellent directionality, which is well-preserved after down-conversion, as the initial and down-converted distributions are comparable. For a heavy mediator, despite the directionality of the initial QPs, the down-converted distribution is much closer to isotropic, particularly for heavier DM. This is simply because heavy mediators and heavy DM favor larger deposits, which lead to a larger number of relaxation events. The directional correlation of the final states with the DM axis is reduced with each relaxation event, so directionality is best preserved when ω/∆ − 2 is small. For this reason, the dashed curves in the bottom-right panel of Fig. 1 show the angular distribution resulting from including only events with total deposits ω < 20∆. These dashed curves retain directionality associated with the low-energy part of the initial angular distribution. Figure 2 shows the estimated sensitivity for the detection of the directional DM wind, with impressive reach. Critically, since the target is isotropic, there is no modulation in the overall rate, unlike many directional detection schemes. Instead, the reach must be defined in terms of the anisotropy in the distribution of the final-state excitations. The dashed lines in Fig. 2 show the projected reach for a detector which counts final-state quasiparticles in each of two angular bins: the 'on-axis' bin, with |cos θ| > 1 2 , and the 'off-axis' bin, with |cos θ| < 1 2 . Such a counting experiment requires a certain number of events in order to establish that the signal is not isotropic at a certain confidence level. 2 In turn, this number of events translates to a minimal cross section for a fixed experimental exposure. Note that for a heavy mediator, we impose a cut to include only small deposits, which reduces the overall rate, but lessens the impact of downconversion.

RESULTS
The two-bin configuration is the minimal experimental configuration for directional detection, and represents the most conservative projection. For heavy mediators, a more ambitious projection is obtained by assuming precise resolution of cos θ, as well as sufficient time resolution so that final-state quasiparticles can be grouped by event.
In this case, the detailed shape of the angular distribution can be compared with the distribution under the null hypothesis of an isotropic background. This again translates to a minimal number of events to detect directionality, and a corresponding minimal cross section, shown by the dotted line in Fig. 2. Further details are given in the SM. For light mediators, since directionality is manifest, this procedure gives almost exactly the same result as the two-bin configuration, so the corresponding line is not shown.
While these are both schematic estimates, the dashed and dotted curves in Fig. 2 are representative of the directional sensitivity of an experiment in the simplest and most sophisticated configurations. Figure 2 shows that a kg-yr exposure of such an experiment would be capable of detecting directionality for DM masses keV m χ GeV, throughout the discovery range of superconducting detectors. Both configurations discussed here would be directionally sensitive at cross sections covering important cosmological targets that are not currently probed by any direct detection experiments.

DISCUSSION
In this work, we have shown that the parameter space to be accessed by superconducting targets [5,6,53] can be probed directionally, even in the case of an isotropic medium. This has important implications for the realworld design of directional superconducting detectors. Directionality will require detectors to push resolutions lower, with thresholds close to the superconducting gap allowing the best coverage of DM parameter space. Detectors that can trap primary QPs are preferred as they will be able to take advantage of the directional correlations in the DM signal. This is in contrast to the weak directionality in the phonon system, in terms of total energy.  [13]. Shaded gray regions indicate existing constraints from SENSEI [50], SuperCDMS HVeV [52], DAMIC [48], Xenon10 [14], DarkSide-50 [44], and Xenon1T [49].
Importantly, in order to scale detectors to kg-yr exposures while retaining directional sensitivity, a massive multiplexing scheme will likely be required. As summarized in Ref. [6], typical detector volumes will be of order 1 cm 3 or smaller in order to attain high collection efficiencies. Mean QP lifetimes are discussed in more detail in Ref. [67], and intrinsic lifetimes (even in the 'dirty' limit) do diverge for very low temperatures and QP occupancy.
For realistic applications of bulk superconducting targets, characterization of mean QP diffusion length in real samples-which will probe the effects of impurities and lattice defects on the scattering lifetime of otherwise ballistic QPs-will determine which materials are best suited for directional DM detection. Such work has been done for large niobium (Nb) crystals [68], but has largely been put aside over the last few decades. These programs will need to be restarted in order to characterize samples with the appropriate properties to detect small energy deposits in the regime relevant to directional DM detection. For materials with a known diffusion length, this directionality can be converted to a rate modulation by making a detector in which one path length to the sensor is much shorter than this diffusion length, and the orthogonal path length is much longer than the diffusion length. Chemical Vapor Deposition-grown superconducting crystals of Nb or Al instrumented on their large surface, with cross-sections of a few mm 2 and thickness of around 100 microns, would achieve this behavior for typical diffusion lengths of a few hundred microns.
On the theoretical side, our results demonstrate directional detection of DM in a target that is otherwise isotropic, in contrast to most directional studies which utilize anisotropy in the material response or in the geometry of the experiment. Here, directionality is inferred from the geometric properties of the excitations themselves, rather than from a rate variation. Our results strongly suggest that for a gapless material with typical acoustic phonon modes, no directional correlation is preserved between the initial DM scatter and the outgoing excitations, as phonons can always be emitted in the limit ∆ → 0. In an anisotropic material, some correlations may still persist further above the gap due to an increased number of forbidden transitions. As an example, indirect-gap materials, if such superconductors exist, would be much more likely to preserve directionality even in the limit of a small gap if large energies are required for inter-valley transitions. We leave the exploration of such scenarios for future work [69], as they require detailed study of the anistropy in the DM response as well.
In this Supplemental Material, we first present the details of quasiparticle (QP) production by dark matter (DM) scattering. We then provide a pedagogical discussion of QP down-conversion in superconductors, supply derivations of the directionality results in the main text, and compare the results to the normal metal case. Finally, we briefly demonstrate the dependence of post-down-conversion directionality on the energy deposited, and show the resulting asymmetry in the recoil spectra.

I. DARK MATTER INTERACTIONS
In this section, we detail our treatment of DM-electron interactions in the language of QP pair production. In particular, we compute Im BCS using the BCS coherence factor F BCS (p 1 , p 2 ), and we show that the free electron scattering picture is recovered in the limit of large deposits. Throughout this section, q denotes the 4-momentum transfer (ω, q).
Consider a DM-electron interaction mediated by a scalar particle ϕ, with interaction Lagrangian where χ is a spin-1/2 DM fermion and ψ is the electron. Using the projection operator P = (1 + γ 0 )/2, we can project out the so-called 'large part' [70] of the electron field ψ s , where s =↑, ↓ refers to the two different spin states. This gives, at lowest order, the interaction Hamiltonian for the electron-mediator interaction: The DM-mediator interaction is governed by a Hamiltonian of the same form with the replacement ψ → χ. Higherorder terms are discussed by Ref. [71]. (Similarly, as shown in that reference, for a light vector mediator A µ , we will have in the low-energy limit L int = g χ A 0 χ † χ + g e A 0 ψ † ψ. Higher-order terms and magnetic interactions will be suppressed by factors of v χ /c in the low-energy limit. In the case of a heavy vector mediator, the A 0 interaction will again dominate, because the currents in the interaction A · j will be suppressed by factors of v χ /c. Thus, light and heavy vector mediators should also be described by interaction Eq. (S.2) in this limit. 3 ) Defining the density operator ρ e (x) ≡ s ψ † s (x)ψ s (x), its Fourier transform is where c p,s annihilates an electron with momentum p and spin s. A similar expression is obtained for ρ χ . We can then write the interaction Hamiltonian as One might object that a propagating A 0 should be suppressed in the non-relativistic limit, by virtue of the constraint ∂µA µ = 0. Indeed, one can compute the time-ordered propagator dt d 3 x exp(iq 0 t − iq · x) 0|T {A µ (x, t), A ν (0, 0)}|0 in the interaction picture and verify that it is not Lorentz-covariant, and that its 00 component is highly suppressed when |q| m A . However, the absence of a kinetic term for A 0 in the Lagrangian introduces an additional terms in the Hamiltonian to precisely cancel this suppression. The effective propagator becomes the Lorentz-covariant propagator −i/(q 2 − m 2 A + iε)(η µν − q µ q ν /m 2 A ), which is unsuppressed for q 0 |q|, or when coupling to a conserved current. This is equivalent to the statement that an off-shell vector can be polarized in any direction. See Sec. 6.2 of Ref. [72] for further insight.
At second order in perturbation theory, the S-matrix will therefore contain a term where ρ e (q) ≡ ρ e (q, ω) = dt e iωt ρ e (q, t) = dt e iωt e iH0t ρ e (q) e −iH0t , (S. 6) with H 0 the free Hamiltonian. In the presence of the lattice potential, an effective electron-electron potential is induced through a phonon loop [73]. The energy eigenstates in the presence of this effective potential are now given by the Bogoliubov QPs, with creation/annihilation operators γ † , γ respectively [58]. To implement the unitary transformation to the QP basis, we simply replace where the coefficients u p and v p satisfy with E and E QP defined as in Eq. (1). We can then isolate the term in ρ e that creates two quasiparticles (breaks a Cooper-pair): giving, according to Eq. (S.6), Plugging this into the S-matrix of Eq. (S.5), and using the fact that the scattering rate is given by Γ = d dt f | f |Ŝ|i | 2 , we find that the lowest order Cooper-pair breaking rate is then simply given by where ω p 1 +p 2 = (p 1 + p 2 ) · v χ − (p 1 + p 2 ) 2 /2m χ is the energy deposited. The quantity is the BCS coherence factor 4 [60]. The E p 1 E p 2 term can be dropped if p 1 or p 2 (or both) are integrated over, provided the remainder of the integrand is even in p 1 and p 2 . But if e.g. p 2 is fixed to be q − p 1 , then this term must be kept when integrating over p 1 . It is straightforward to see that this matches onto the rate for free electron scattering when p 1 and p 2 are away from the gap (|p 2 i /2m * − E F | ∆). In this limit, because of the functional form of the coherence factor, we always have That is, the kinematical constraint reduces to that of ordinary non-relativistic scattering. Meanwhile, the quantity is known in the literature as the dynamic structure factor, in terms of which we have, at lowest order in perturbation theory, In our case, the structure factor is simply given by allowing the loss function, Im(−1/ BCS ), to be expressed in terms of the QP dispersion relation and the BCS coherence factor. When higher order terms are included, we can resum the series (at zero temperature) in the random phase approximation (RPA) [74,75]:Ŝ where (RPA) (q, ω) ≡ 1 + χ(q, ω), for 1 − χ(q, ω) ≡ −1/ (1) (q, ω). Importantly, we have Im (RPA) (q, ω) = Im −1 (1) (q, ω) = Im χ(q, ω) = πe 2 q 2 S(q, ω) .

(S.18)
Putting everything together, we have 19) and writing the coherence factor in terms of (1) gives . (S.20) To the accuracy that (RPA) represents the true dielectric function , we have derived Eq. (2). 5 The full form of (RPA) in the BCS vacuum will be discussed in future work [61]. In this Letter, we make the approximation where the Lindhard function [59,74] FEG is the RPA dielectric function for a free electron gas (FEG), and accounts for screening and in-medium effects in a normal metal. Fortunately, by Eq. (S.18), Im This is because the imaginary part depends only on the spectrum of the Hamiltonian and on its relation to the operator ρ e . It corresponds precisely to the sum over states in the scattering rate.

II. QUASIPARTICLE DOWNCONVERSION
We simulate QP down-conversion following a similar procedure to the calculations of Refs. [65,67,77,78]. The principal difference is that we track the full momentum vector for scattered particles to retain information about the scattering direction relative to the momentum of the initial scattering event. Here we review the relevant scattering rate calculations, and we derive the emission angles of QPs and phonons in each relaxation process. We compare the final result to the well-established normal metal case described in Refs. [79,80]. Down-conversion in the limit of high-energy initial QPs is also discussed by Refs. [65,66].

A. Phonon Scattering at the Fermi Surface
We treat QP-phonon interactions following Ref. [67]. The simplified model in this treatment contains a single acoustic phonon branch with the dispersion relation ω q = qc s , where q is the phonon momentum and c s is the sound speed in the material. For the dynamics of the problem, scattering is contained to the first Brillouin zone (q < 2π a ), so we do not include an explicit upper limit in momentum in the rate integral. As significant down-conversion already occurs at energies well below the optical phonon modes in most superconductors, we do not explicitly include optical phonon emission in our calculations. This is compatible with the conclusions of past down-conversion codes (see e.g. Refs. [65,77]).
For acoustic phonon emission, we first adopt the result of Ref. [67] for the emission rate of an acoustic phonon of momentum q in the zero-temperature limit: where Z 0 is the 'renormalization parameter' defined in Ref. [67] (Typically Z 0 ∼ 2). We have also assumed that α 2 (q)F (q) ≈ b ω 2 q , the Debye solution for this quantity [79]. The QP energy E QP is defined in Eq. (1).

B. Computing the Scattering Angle
Given the general scattering rate of Eq. (S.23), we also need to determine the scattering angle for a phonon of energy ω q = qc s emitted by an QP of energy E QP (k). We find this angle by first solving for k from the conservation of energy relation, E QP (k ) = E QP (k) − qc s , where E QP (k) ≡ E QP (k) as defined in Eq. (1). Solving, we get In a normal metal or a semiconductor, the p F dependence cancels. However, for a superconductor with ∆ > 0, the p F dependence is retained.
We can now use momentum conservation to solve for scattering angle. Writing k = k + q, we have k 2 = k 2 + q 2 − 2kq cos θ q . Solving explicitly for k in the QP dispersion relation gives where s = ±1. An analogous sign s appears in the solution for k . We then find We now consider a few limiting cases to understand the angular spread in the phonon spectrum. First, observe that for QPs far from the gap, with ∆ E QP E F , the scattering angle is unrestricted. In this limit, we have For typical materials, γ is O(1) and E F ∆. Thus, if s = s , the last term dominates in the limit of small ω q , and | cos θ q | = 1 is allowed. On the other hand, if s = s , then the last term vanishes. The first term can be made arbitrarily small in the small-ω q limit, and can even be made to cancel with the second term, which is itself always small. In this case, cos θ q = 0 is allowed.
On the other hand, consider the near-gap regime, where k ∼ p F . Here we can write E QP = (1 + δ)∆ with δ 1, and since the final-state QP energy is at least ∆, we must have ω q = aδ∆ with 0 < a < 1. Inserting these expressions into Eq. (S.26) and expanding for small δ gives (S.28) However, the phonon emission process is kinematically forbidden if |cos θ q | > 1, and minimizing Eq. (S.28) over a gives min a∈(0,1) Thus, phonon emission is only kinematically allowed for sufficiently large δ, i.e., for This result is self-consistent: in typical materials, γ is O(1), and E F ∆, so δ min 1. This gives rise to a condition for phonon emission: QPs with energies below this threshold are ballistic: no phonon emission is allowed. The angular distribution of final-state phonons is peaked oppositely to that of final-state QPs, and is closer to the isotropic distribution. The distribution is shown explicitly in Fig. S2 for the same cases as in Fig. 1.

C. Relation to Normal Metal Scattering
It is instructive to compare Eq. (S.23) to the equivalent rate in the normal metal case. In a normal metal, this emission rate becomes [79] Γ(E k ) = 2π where α q is the coupling for electron-phonon scattering. If we assume that scattering is isotropic, we can make a change of variables such that dp ∝ dω q d(cos θ), where cos θ is the scattering angle between k and k . For scattering near the Fermi energy in the metal, conservation of momentum gives where p F = √ 2mE F is the Fermi momentum. This implies d(cos θ) = −q dq/p 2 F , which in turn allows us to write the differential scattering rate for ω q < E QP as Here the subscript 'n' indicates the normal metal case; N 0 is the normal metal density of states at E F ; and α 2 F (q) is the coupling-weighted phonon density of states, given by For the scaling used earlier this gives the normal result that the acoustic scattering rate goes as E 3 QP when integrated over energy, and the differential spectrum goes as ω 2 q . When we derive the emission rate for superconductors, there are two important modifications required to get from Eq. (S.32) to the final differential rate in emitted phonon energy. First, we modify the dispersion relation to that of the QPs in the superconductor. In the metal, we had E k = E k ≡ p 2 /(2m * ) − E F relative to the Fermi surface, but in a superconductor, the gap energy modifies this to (S. 36) In the normal metal phonon interaction, the total coupling has a density of states term that is valid in the metal, but not in the superconductor, since there are no states at E k < ∆. We can use the fact that the total number of states is the same to find the modified density of states at a given energy, i.e., we have dE N s (E) = dE N n (E), where the subscript 's' indicates the superconductor case. We thus find that Thus, to rescale α 2 F (q) for the superconducting case, we have to rescale by the ratio of superconducting to normal states at a given energy. This gives us the first additional factor in the superconducting rate equation. Ref. [60] points out that, in principle, the gap function is complex-valued, hence the need to take only the real part.
The second correction factor is the coherence factor described in the previous section. This is a purely BCS effect that arises from the collective nature of the superconducting states. This factor ensures that the divergence in the density of states does not lead to a divergence in the phonon scattering rate. Taking the type I coherence factor for phonon emission from Ref. [60], we find using the same normalization procedure as Ref. [67]. Combining this with the previous correction factor and multiplying by the normal metal scattering rate produces the scattering rate of Eq. (S.23). This heuristic argument elucidates the origin of Eq. (S.23) in the low-temperature limit with a less formal approach than that in Ref. [67].

III. ESTIMATION OF DIRECTIONAL REACH
In this section, we detail the methodology used to produce Fig. 2 in the main text, and we demonstrate the impact of total-deposit cuts on the directionality of the signal.

A. Statistical methods
The directional reach in Fig. 2 is estimated in two distinct ways. The dashed curves are based on measurement of an asymmetry in the counts of final-state QPs between two bins of equal solid angle, and the dotted curve is based on comparison of the full angular distribution against the null hypothesis.

Two-bin reach
We first discuss the two-bin reach estimate. The premise of this test is that an isotropic background gives rise to an isotropic distribution of final-state quasiparticles, and, in particular, produces statistically-indistinguishable counts in any two bins of equal solid angle. Since the DM wind is not isotropic, it is possible to statistically distinguish the counts in the two bins given a sufficient number of events.
We make the simplifying assumption that the angles of the final-state QPs are independent random variables. This is not strictly the case, since QPs that originate from the same event have some angular correlation. However, given a large number of events, such correlations are extremely sparse. Moreover, by simulating an ensemble of isotropic DM scattering events, we have directly checked that such correlations are irrelevant at the number of events needed to establish directionality. Having made this assumption, the assignment of an angular bin to each QP can be treated as a Bernoulli trial. We can then use the binomial test to determine whether to reject the isotropic distribution given a particular set of QPs.
This procedure allows us to determine whether a particular sample of final-state QPs is consistent with an isotropic signal. Next, we must translate this to a minimum number of events needed to establish directionality. To that end, we randomly draw samples of N QP = 2, 3, 4, . . . QPs and evaluate the binomial test for each sample, repeating the process many times for each fixed N to obtain a median p-value. We advance N QP until this median p-value drops below the threshold value of 0.05, and we interpret the resulting value of N as the typical number of QPs needed in order to establish directionality. Finally, this number of QPs must be translated to a number of scattering events. We estimate this number as N e ≡ N QP /n, wheren is the average number of final-state QPs produced by a scattering event. This N e is indicated by the dashed curves in Fig. 2.
For the heavy mediator case, we introduce an additional step. As we discuss below, the distribution of finalstate QPs produced by an event approaches the isotropic distribution as the deposit becomes much larger than the superconducting gap 2∆. Thus, it is advantageous to restrict attention to events with total deposit below some cut, even at the cost of a reduced event rate below the cut. The dashed curve in the right panel of Fig. 2 is a composite of two reach curves obtained with cuts ω < 10∆ and ω < 50∆. Due to the complicated relationship between the total deposit, the initial directionality, and the effects of down-conversion, each of these cuts preserves overall directionality for a different mass range, and the combination of the two gives directional sensitivity over the entire mass range. Due to the large deposits favored in the heavy-mediator case, a cut on the total deposit is essential to establish directionality for all but the lowest masses.
This simplistic treatment produces a rough upper bound on the number of events needed to detect directionality, and admits a very direct interpretation. Even at the level of a two-bin experimental configuration, more sophisticated statistical treatments may yield slightly stronger results. In particular, it is possible to extract directionality using the Skellam distribution, as in rate modulation experiments [38]. Here one treats the count in each bin as a Poisson random variable, so that the difference in the number of counts between the two bins is a Skellam-distributed random variable. It is then possible to test whether the two Poisson counts are produced with the same rate. However, our case is simpler than a traditional rate modulation experiment in that the target is isotropic, so the total rate is fixed. We have checked that using the Skellam distribution offers a slight enhancement to the two-bin reach, but does not qualitatively change the result. Figure 2 also includes an estimate based on the full angular distribution of the final-state QPs, assuming an experimental configuration with great angular precision. For this estimate, we begin with an ensemble of simulated DM scattering events. Next, a second 'null' ensemble of DM scattering events is simulated with an isotropic distribution of DM directions, using the same speed distribution as the Standard Halo Model. For each simulated event, we compute the mean angle of the final-state QPs, cos θ i . We thus obtain two sets of samples { cos θ i } SHM and { cos θ i } null for the two ensembles. We then determine the average number of events needed to reject at 95% C.L. the hypothesis that Left: since ω = 2.09∆, neither initial QP can emit a phonon with energy above 2∆, so no phonons can produce QP pairs. Thus, the only QPs in the final state are those produced in the DM scattering event, with their directions barely altered. Right: now ω = 12.44∆ and the emission of above-gap phonons is allowed. Thus, a chain of phonon emissions and decays erases much of the initial directionality, although a preference for off-axis final states remains visible.

Full angular distribution
the 'SHM' and 'null' samples are drawn from the same distribution, using the two-sample Kolmogorov-Smirnov test. Note that unlike the binomial test of the previous case, the samples being compared in this case are truly independent: since each mean angle corresponds to a single event, and vice versa, all of the values we draw originate from different scattering events, and they are thus independent random variates. For light mediators, the two procedures give a nearly identical result, and Fig. 2 shows only the simpler two-bin result. This is easily understood in light of the typical number of QPs produced in a scattering event: for a light mediator, a typical final state consists of an O(1) number of QPs, and deviation of these QPs from the isotropic distribution generally gives rise to a two-bin asymmetry. In fact, for large DM masses, taking the mean angle in each event discards information to the point that the distributional test slightly underperforms the two-bin test. For a heavy mediator, on the other hand, a typical final state may consist of O(10 4 ) QPs for large DM masses. The angular mean cos θ i will typically be very close to zero, and distinguishing the distribution of these means from the isotropic case becomes a problem of precision measurement. The large number of QPs makes this a realistic possibility. However, we estimate that achieving the high-resolution dotted curve in Fig. 2 would require a measurement precision of O(10 −2 ) in cos θ i .

B. Energy dependence of directionality
In this section, we demonstrate the dependence of final-state directionality on the deposited energy. As noted in the main text, larger deposits allow for a larger number of relaxation events during the down-conversion process, which attenuates the correlation between the directions of the final-state excitations and those of the initial excitations produced by DM scattering. This is illustrated for particular realizations in Fig. S3. In the left panel, the relatively small number of relaxation events and the low energies of the emitted phonons guarantee that the directions of the initial QPs are well-preserved, and no additional QPs are produced. In the right panel, on the other hand, the Left: asymmetry A2 of final-state QPs (blue) and phonons (orange). Directionality of the final excitations is quickly erased for ω ∆, and A2 approaches zero in the large-ω limit. Right: Fraction of the total energy residing in the QP system after down-conversion. The gray band shows the range of asymptotic results 0.57-0.60 obtained in the literature for bulk Al superconductors [65,77].
larger deposit allows for a larger number of relaxation events, with additional QP pair production. While there remains a directional correlation between the final state and the initial QPs, this correlation is partially erased by down-conversion. Directional information is effectively lost for very large deposits.
The impact of down-conversion means that small deposits are favorable for directionality even if the directionality of the initial QPs is smaller in this regime. We can demonstrate this explicitly by evaluating the final-state asymmetry as a function of the deposited energy. To facilitate quantitative discussion of directionality, we introduce a quantitative 'two-bin asymmetry' A 2 , defined as follows. As discussed above, we divide the final-state QPs into two bins of equal solid angle: the 'on-axis' bin, with |cos θ| > 1 2 , and the 'off-axis' bin, with |cos θ| < 1 2 . We denote the counts in these two bins by n on and n off , respectively, and then define A 2 ≡ 2 × n on n on + n off − 1 .
(S. 39) In particular, the isotropic distribution gives A 2 = 0, and a totally asymmetric distribution gives A 2 = 1.
The left panel of Fig. S4 shows A 2 for the final-state QPs produced by a single QP injected at fixed energy ω with cos θ = 1. This serves as a proxy for the preservation of directionality at fixed deposit. Directionality is almost perfectly preserved for very low-energy QPs with ω < 3∆. Above this threshold, it becomes possible for the QP to emit an above-gap phonon, with E ph > 2∆. Such a phonon subsequently decays to another pair of QPs, which have only weak angular correlation with the original QP. The nature of the 3∆ threshold is also clearly visible in the right panel of Fig. S4, which shows the fraction of the total energy that resides in the QP system after down-conversion. For E QP < 3∆, the QP relaxes almost all the way to the gap by emission of sub-gap phonons, which cannot produce any additional QPs. Thus, the final state consists of a single QP with E QP ≈ ∆, and a set of phonons with all the remaining energy from the deposit. The fraction of the initial energy in the QP system is approximately ∆/ω. Upon reaching ω > 3∆, emission of above-gap phonons produces additional QPs in the final state, sharply raising the fraction of the total energy in the QP system. At large ω, this fraction asymptotically reaches ∼ 0.6. This is consistent with previous studies of down-conversion in the high-energy limit, which find fractions between 0.57 and 0.60 [65,77] (gray band in Fig. S4).
Naively, Fig. S4 suggests that the smallest deposits will yield the strongest directionality. However, imposing an upper limit on the deposit also influences the directionality of the initial excitations, and, of course, the rate of events which fall below the cut. To study the total directionality as a function of the deposited energy, in Fig. S5, we show the spectrum of final-state QPs in each of the two angular bins ('on-axis' and 'off-axis') for several DM masses in the light-and heavy-mediator limits. The bottom panels of Fig. S5 show the ratios of these spectra, i.e., a signed and Top panels: total event rates as a function of initial deposit in the on-axis (|cos θ| > 1 2 , solid) and off-axis (|cos θ| > 1 2 , dashed) bins for several DM masses, for light (left panel ) and heavy (right panel ) mediators. Normalization of each pair of curves is arbitrary and fixed for ease of visualization. Bottom panels: ratio of off-axis to on-axis QP counts. The left and right panels assume a light and heavy mediator, respectively.
shifted version of the two-bin asymmetry A 2 .
For light mediators, directionality is quickly lost for deposits well above the gap, and the ratio approaches 1. Moreover, the ratio is generally above 1. For heavy mediators, due to the directionality of the initial excitations, the ratio declines less noticeably, but it is near 1 throughout the plot and asymptotically reaches 1. As anticipated in Fig. 1, light mediators always enhance the off-axis rate, while heavy mediators can enhance either the off-axis or on-axis rates, depending on the DM mass: since scattering through a heavy mediator can produce a QP distribution peaked either in the forward direction or off-axis, the ratio can be either below or above 1.