Directional Detectability of Dark Matter With Single Phonon Excitations: Target Comparison
Ahmet Coskuner, Tanner Trickle, Zhengkang Zhang, Kathryn M. Zurek
CCALT-TH-2021-008
Directional Detectability of Dark Matter With Single Phonon Excitations:Target Comparison
Ahmet Coskuner,
1, 2
Tanner Trickle, Zhengkang Zhang, and Kathryn M. Zurek Department of Physics, University of California, Berkeley, CA 94720, USA Theoretical Physics Group, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA Walter Burke Institute for Theoretical Physics, California Institute of Technology, Pasadena, CA 91125, USA
Single phonon excitations are sensitive probes of light dark matter in the keV-GeV mass window.For anisotropic target materials, the signal depends on the direction of the incoming dark matterwind and exhibits a daily modulation. We discuss in detail the various sources of anisotropy, andcarry out a comparative study of 26 crystal targets, focused on sub-MeV dark matter benchmarks.We compute the modulation reach for the most promising targets, corresponding to the cross sectionwhere the daily modulation can be observed for a given exposure, which allows us to combine thestrength of DM-phonon couplings and the amplitude of daily modulation. We highlight Al O (sapphire), CaWO and h-BN (hexagonal boron nitride) as the best polar materials for recoveringa daily modulation signal, which feature O (1 – 100)% variations of detection rates throughout theday, depending on the dark matter mass and interaction. The directional nature of single phononexcitations offers a useful handle to mitigate backgrounds, which is crucial for fully realizing thediscovery potential of near future experiments. I. INTRODUCTION
If the cold dark matter (DM) in the universe consistsof new particles, they must interact very weakly withthe Standard Model. Directly detecting these feeble in-teractions in a laboratory requires extraordinarily sen-sitive devices. Traditional direct detection experiments( e.g.
ANAIS [1], CRESST [2–4], DAMA/LIBRA [5],DAMIC [6], DarkSide-50 [7], DM-Ice [8], KIMS [9],LUX [10, 11], SABRE [12], SuperCDMS [13, 14], andXenon1T [15]), based on nuclear recoil, are graduallyimproving their sensitivity and closing the open param-eter space before reaching the irreducible solar and at-mospheric neutrino background. However, these ex-periments are fundamentally limited in the DM mass, m χ , they can probe. When the DM scatters off a nu-cleus at rest, the energy deposited, ω , is limited by ω (cid:46) m χ v /m N , with v ∼ − , and vanishes quicklyas the DM mass decreases below the nucleus mass m N .This limitation in DM mass is typically not problem-atic in the search for the prototypical weakly interact-ing massive particle (WIMP) which produces the DMabundance through freeze-out, as m χ (cid:46) GeV would bothbe overabundant and be in tension with indirect detec-tion bounds on DM annihilation rates. However, manyother theoretically motivated explanations of the originof DM such as freeze-in [16, 17], hidden sector DM [18–21], asymmetric DM [22–25], and strong self interactions[26, 27], allow for DM ligher than a GeV and thereforeshould be searched for by means other than nuclear re-coil.In the pursuit of sub-GeV DM, several new experimen-tal concepts have been proposed. These include electronexcitations in a variety of target systems [28–40] for DMwith mass above an MeV, while phonon excitations [41–55], with energies up to O (100) meV, have been shown tobe especially sensitive to a wide range of DM models with masses down to a keV. This coupled with the fact thatdetector energy thresholds are approaching the O (100)meV range [56–60] makes phonon excitations an excitingavenue for DM direct detection. Phonons are quasiparti-cle vibration quanta which can exist in multiple states ofmatter, e.g. as sound waves in liquids or superfluids andlattice vibrations in crystalline solids. Superfluid heliumhas been proposed [41] and studied as a light DM detec-tor [42–46] and an experiment is currently in the R&Dphase [61]. Crystal targets have also been proposed [47]and studied extensively. Initial studies focused on GaAsand Al O (sapphire) targets [48], and more recently thisanalysis has been extended to account for more generalDM interactions [49, 53] and applied to a broader setof target materials [50]. Other targets have also beenproposed individually, e.g. SiC [52], and there has beenwork on understanding the signal from multi-phonon ex-citations [51, 54, 55]. Similar to superfluid helium, a DMdetector using a crystal target with phonon readout isalso in the R&D phase of development [61].Most of the previous work has focused on calculat-ing the theoretically predicted DM-phonon interactionstrength. Equally important is to minimize the experi-mental background. This becomes easier when the DMscattering signal has unique properties which can distin-guish it from backgrounds. For example, in experimentssensitive to nuclear recoil or electron excitations, the ratemodulates annually due to the change in the DM veloc-ity distribution in the Earth frame, as the Earth orbitsaround the Sun [30, 31, 62].In a phonon-based experiment, the DM scattering ratecan have a larger and more unique signature: daily mod-ulation. As the Earth rotates about its own axis, the ori-entation of the detector relative to the DM wind changes.In a nuclear recoil experiment this does not have an ef-fect since the interaction matrix element is independentof the direction of the DM velocity relative to the detec-tor orientation — an (unpolarized) nucleus is isotropic a r X i v : . [ h e p - ph ] F e b in its response. However, crystal targets can be highlyanisotropic, which means that the amplitude of the re-sponse depends not only on the magnitude of the mo-mentum transfer but also on its direction. This can leadto a significant daily modulation. Moreover, since themodulation pattern depends on the crystal orientation,running an experiment with multiple detectors simulta-neously with different orientations can further enhancethe signal-to-noise ratio. This effect was studied for sap-phire in Ref. [48] (see also Refs. [33, 49, 63, 64] for discus-sions of daily modulation in electron excitations). In thiswork, we expand the understanding of the daily modula-tion effect in single phonon excitations to a broader rangeof materials, including those targeted in Ref. [50].In particular, we highlight the following targets in themain text: Al O (sapphire) and CaWO , already uti-lized for DM detection and have a significant daily mod-ulation; SiO , shown to have a strong reach to severalbenchmark models; SiC, proposed in Ref. [52]; and h-BN (hexagonal boron nitride), a highly anisotropic ma-terial. Among them, Al O and CaWO have the bestprospects overall in terms of daily modulation reach andexperimental feasibility. Of the additional materials con-sidered in Ref. [50], results for those with daily modula-tion larger than 1% (for at least some DM masses) arepresented in an appendix. We make available our codefor the phonon rate calculation [65], and also publish aninteractive webpage [66] where results for all the materi-als presented Ref. [50] and in this paper, including reachcurves, differential rates and daily modulation patterns,can be generated from our calculations. II. DIRECTIONAL DETECTION WITHSINGLE PHONON EXCITATIONSA. Excitation Rate
We begin by summarizing the formulae for singlephonon excitation rates; see Refs. [47–49] for more de-tails. For the scattering of a DM particle χ with mass m χ and general spin-independent interactions, the rateper unit target mass takes the form R ( t ) = 1 ρ T ρ χ m χ πσ ψ µ χψ (cid:90) d v f χ ( v , t ) × (cid:90) d q (2 π ) F ( q ) S (cid:0) q , ω q ) , (1)where v is the incoming DM’s velocity, q is the momen-tum transferred to the target, ρ T is the target’s massdensity, and ρ χ = 0 . is the local DM density. Twenty-six materials are included on the interactive web-page [66]: Al O , AlN, h-BN, CaF , CaWO , CsI, C (diamond),GaAs, Ge, GaN, GaSb, InSb, LiF, MgF , MgO, NaCl, NaF, NaI,PbS, PbSe, PbTe, Si, SiC, SiO , ZnO, and ZnS. σ ψ , with ψ = n or e (neutron or electron), is a referencecross section defined as σ ψ ≡ µ χψ π |M χψ ( q = q ) | , (2)where µ χψ is the reduced mass, M χψ is the vacuum ma-trix element for χψ → χψ scattering, and q is a referencemomentum transfer. We present the reach in terms of σ ψ ,with q = m χ v (where v = 230 km/s, the dispersion ofDM’s velocity distribution) for ψ = n and q = αm e for ψ = e . f χ ( v , t ) is the DM’s velocity distribution in thelab frame, taken to be a truncated Maxwell-Boltzmanndistribution, boosted by the time-dependent Earth veloc-ity v e ( t ), as will be discussed in more detail in the nextsubsection. F med ( q ) is the mediator form factor, whichcaptures the q dependence of the mediator propagator: F med ( q ) = (cid:40) , ( q /q ) (light mediator) . (3)Finally, S ( q , ω ) is the dynamic structure factor that en-codes target response to DM scattering with momentumtransfer q and energy transfer ω , constrained by energy-momentum conservation to be ω q = q · v − q m χ . (4)Generally, one sums over a set of final states f with en-ergies ω f , and S ( q , ω ) takes the form S ( q , ω ) = (cid:88) f π δ ( ω − ω f ) S (cid:48) f ( q ) . (5)For single phonon excitations, we assume the target sys-tem is initially prepared in the ground state at zero tem-perature with no phonons, and sum over single phononstates labeled by branch ν and momentum k inside thefirst Brillouin zone (1BZ). Lattice momentum conserva-tion dictates that q = k + G , with G a reciprocal latticevector. To find k and G from a given q , we first find thereduced coordinates ( q , q , q ) ( i.e. q = (cid:80) i =1 q i b i with b i the basis vectors of the reciprocal lattice), and thenfind the nearest point ( G , G , G ) with G i ∈ Z . In thisway, any q outside of the 1BZ is mapped to a k inside the1BZ and a G vector. The sum over final states thereforeonly runs over ν , S ( q , ω ) = (cid:88) ν πδ ( ω − ω ν, k ) S (cid:48) ν ( q ) . (6)As was shown in Refs. [47–49], S (cid:48) ν can be written in termsof the phonon energies ω ν k , eigenvectors (cid:15) ν k j and an ef-fective DM-ion couplings Y j (with j labeling the ions inthe primitive cell): S (cid:48) ν ( q ) = 12Ω ω ν, k (cid:12)(cid:12)(cid:12)(cid:12)(cid:88) j e − W j ( q ) √ m j e i G · x j (cid:0) Y j · (cid:15) ∗ ν, k ,j (cid:1)(cid:12)(cid:12)(cid:12)(cid:12) , (7)where Ω is the volume of the primitive cell, and m j , x j ,and W j ( q ) ≡ Ω4 m j (cid:80) ν (cid:82) d k (2 π ) | q · (cid:15) ν, k ,j | ω ν, k are the masses,equilibrium positions, and Debye-Waller factors of theions, respectively. We obtain the material-specific forceconstants in the quadratic crystal potential and the equi-librium positions from density functional theory (DFT)calculations [48, 50, 67], and use the open-source phononeigensystem solver phonopy [68] to derive the values of ω ν, k , (cid:15) ν, k ,j for each material.The DM-ion coupling vectors Y j are DM model de-pendent. In our target comparison study in Sec. III, wewill focus on two sets of benchmark models, with a lightdark photon mediator and a heavy or light hadrophilicscalar mediator, respectively. These are the same modelsconsidered in Ref. [50], for which Y j are given by Y j = − q · Z (cid:63)j ˆ q · ε ∞ · ˆ q (dark photon med.) , q A j F N j ( q ) (hadrophilic scalar med.) . (8)Here Z (cid:63)j is the Born effective charge tensor of the j th ion, ε ∞ is the high-frequency dielectric tensor that capturesthe electronic contribution to in-medium screening, A j isthe atomic mass number, and F N j ( q ) = j ( qr j ) qr j e − ( qs ) / (with r j = 1 . A / j fm, s = 0 . Y j scales with q and the energy conserving deltafunction contributes a factor of q − (see Eq. (11) below),we see that for a heavy mediator, the integral scales as (cid:82) dq q ω − and so is always dominated by large q . For alight mediator, on the other hand, the integral scales as (cid:82) dq q − ω − . So for optical phonons with ω ∼ q , it re-ceives similar contributions from all q , whereas for acous-tic phonons, it is dominated by small q where ω ∼ q . B. Daily Modulation
In the rate formula Eq. (1), the time dependence comesfrom the DM’s velocity distribution f χ ( v , t ), specificallyvia the Earth’s velocity v e ( t ) that boosts the distribution.Concretely, we take f χ ( v , t ) = N exp (cid:104) − ( v + v e ( t )) v (cid:105) Θ (cid:0) v esc −| v + v e ( t ) | (cid:1) , (9)where N is a normalization constant such that (cid:82) d v f χ ( v ) = 1, v esc = 600 km/s is the galactic escapevelocity, and v = 230 km/s as mentioned above. Assum-ing the detector is fixed on the Earth, in the lab frame v e becomes a function of time that is approximately peri-odic over a sidereal day as a result of the Earth’s rotation.As a default setup, we adopt the detector orientation inRefs. [48, 49, 63], for which, independent of the detector’slocation, v e ( t ) = v e sin θ e sin φ ( t )sin θ e cos θ e (cos φ ( t ) − θ e + sin θ e cos φ ( t ) , (10)where v e = 240 km/s, θ = 42 ◦ , and φ ( t ) = 2 π (cid:0) t
24 hr (cid:1) . Itis this periodicity of the direction of v e ( t ) that inducesthe daily modulation in R ( t ) we study in this work. With the specific form of f χ in Eq. (9), the velocityintegral in Eq. (1) can be done analytically [48, 49, 63].We define g ( q , ω, t ) ≡ (cid:90) d v f χ ( v , t ) 2 πδ ( ω − ω q )= π v N q (cid:110) exp (cid:104) − ( v − ( q ,ω,t )) v (cid:105) − exp (cid:104) − v v (cid:105)(cid:111) , (11)where v − ( q , ω, t ) = min (cid:16)(cid:12)(cid:12)(cid:12) ˆ q · v e ( t ) + q m χ + ωq (cid:12)(cid:12)(cid:12) , v esc (cid:17) . (12)We will refer to g ( q , ω, t ) as the kinematic function . Therate formula Eq. (1) then becomes R ( t ) = 1 ρ T ρ χ m χ πσ ψ µ χψ (cid:90) d q (2 π ) F ( q ) (cid:88) ν S (cid:48) ν (cid:0) q ) g ( q , ω ν, k , t ) . (13)With the rate written in this form, the time dependencenow comes from the v − function contained in g ( q , ω, t ).As we will discuss in detail in the rest of this subsection,the origin of daily modulation is as follows. First, the Annual modulation is also present, due to the change of the magnitude of v e , as in any terrestrial direct detection experiment.Here we fix v e = 240 km/s and focus on the daily modulationsignal, which is unique to anisotropic (crystal) targets. - - - - - - - - - - - - - - - - - - - - - - Figure 1.
Top:
To understand the kinematic function, g ( q , ω ), defined in Eq. (11), we plot v ∗ ≡ q m χ + ωq as a function of q (blue) for various m χ and ω values. Comparing v ∗ to v e and v e + v esc we can qualitatively reconstruct the shape of g ( q , ω ),as discussed in the text. Bottom: g ( q , ω ) vs. q for several fixed m χ , ω values, with varying ˆ q · ˆ v e . The kinematic functionweights different ˆ q directions according to their angle with respect to v e ( t ), which ultimately leads to a daily modulating rate. kinematic function g ( q , ω, t ) selects a region of q space ateach time of the day that is strongly correlated with v e ( t ).For anisotropic targets, this then results in a modulatingrate after the q integral in Eq. (13). Intuitively, the DMwind hits the target from different directions throughoutthe day, some of which may induce a stronger responsethan others.
1. Kinematic Function
The kinematic function g ( q , ω ν, k , t ) can be viewed as aweight function: for each phonon branch ν , the integrandin Eq. (13), F ( q ) S (cid:48) ν (cid:0) q ), is weighted toward momen-tum transfers q that maximize the g function or, equiva-lently, minimize v − defined in Eq. (12). To visualize thisminimization, we plot v ∗ ≡ q m χ + ωq (14)as a function of q in the top panel of Fig. 1. Setting ω to a constant approximates the case of optical phonons,which have relatively flat dispersions, whereas ω → ω/q is bounded by the sound speed, which is typicallymuch smaller than the DM’s velocity. We can identify three distinct regions (as shown with different colors inthe plot): • For v ∗ ≥ v esc + v e , we have v − = v esc and therefore g = 0 for all ˆ q directions. This is the kinematicallyforbidden region. • For v e ≤ v ∗ < v esc + v e , the g function is maximizedat ˆ q · ˆ v e = − • For v ∗ ≤ v e , the g function is nonzero for all ˆ q directions, and is maximized at ˆ q · ˆ v e = − v ∗ /v e . Inthe large m χ , small ω limit, v ∗ →
0, and thereforethe g function is maximized when ˆ q · ˆ v e = 0.These behaviors are seen in the lower panels of Fig. 1(see also Ref. [63]), where we plot g ( q , ω ) as a functionof q for fixed m χ , ω , and with varying ˆ q · ˆ v e . Note thatin the ω → g function has support down to q = 0, but the phase space integral for acoustic phononsis cut off at q min (cid:39) ω min c s = 2 × − keV (cid:0) ω min (cid:1)(cid:0) × − c s (cid:1) ,where ω min is the detector’s energy threshold, and c s isthe sound speed (slope of the linear dispersion). Fromthese plots we see that the kinematically favored region of q is strongly correlated with ˆ v e ( t ) and, therefore, rotateswith it throughout the day. This rotation then translatesany target anisotropy into a detection rate that modu-lates daily. Figure 2. Comparison between the various sources of anisotropy a in SiO target, for an example DM mass for each benchmarkmodel. A 1 meV energy threshold is assumed in all cases. As discussed in the text, anisotropy in the Y j · (cid:15) ν, k ,j factor in Eq. (7)is the dominant factor in determining the daily modulation pattern.
2. Sources of Anisotropy
There are a number of possible sources of anisotropy,as we can infer from Eq. (13) and Eq. (7). First of all, thephonon energies ω ν, k generically depend on the directionof q = k + G . This means that the region selected by thekinematic function, as discussed above, does not preserveits shape as it rotates in q space. Also, the ω − ν, k factorin Eq. (7) is different in the dominating kinematic regionat different times of the day, which adds to the dailymodulation signal.The anisotropy in ω ν, k has two contributing factors.First, phonon dispersions can be anisotropic as a re-sult of crystal structures. For example, in h-BN, thesound speed of the longitudinal acoustic phonons dif-fers by more than a factor of two between different k directions. Second, by the prescription explained aboveEq. (6), a sphere of constant q outside the 1BZ does notmap to a sphere of constant k inside the 1BZ. Since thesize of the 1BZ is typically O (keV), and the DM velocityis O (10 − ), this is relevant for m χ (cid:38) MeV. Another re-lated source of anisotropy is the e i G · x j factor in Eq. (7):a constant- q sphere outside the 1BZ does not map ontoa unique G vector.In addition to ω ν, k and e i G · x j discussed above, thescalar product of the DM-ion coupling and phonon eigen-vectors, Y j · (cid:15) ν, k ,j , can also be anisotropic for a vari-ety of reasons, depending on the DM model. For thehadrophilic scalar mediator model, Y j · (cid:15) ν, k ,j are simplyproportional to the longitudinal components of phononeigenvectors ˆ q · (cid:15) ν, k ,j , so the anisotropy is determinedby the extent to which the phonon eigenvectors deviatefrom transverse and longitudinal in different ˆ q directions.For the dark photon mediator model, Y j · (cid:15) ν, k ,j are in- stead proportional to ˆ q · Z (cid:63)j · (cid:15) ν, k ,j ˆ q · ε ∞ · ˆ q , so there are additionalanisotropies if the Born effective charges Z (cid:63)j and dielec-tric tensor ε ∞ are not proportional to the identity. Allthese anisotropies are ultimately determined by the crys-tal structure.We can carry out a simple exercise to see how thevarious sources of anisotropy discussed above contributeto the full daily modulation signal. As an example, weconsider a SiO target, and pick one m χ value for eachbenchmark model, as shown in the three panels of Fig. 2.We obtain the full rate normalized to its daily average, R/ (cid:104) R (cid:105) , as a function of time, as shown by the solid redcurves labeled by “full.” We then artificially make thevarious factors in the rate formula isotropic and see howthe modulation pattern changes.First, we make S (cid:48) ν ( q ) isotropic by setting ω ν, k and Y j · (cid:15) ν, k ,j to their values at a specific direction ( ˆ q = ˆ z ),and setting e i G · x j →
1. This isolates the effect of thekinematic function g ( q , ω ν, k , t ) on daily modulation. Theresults are shown by the dot-dashed purple curves inFig. 2, labeled “isotropic S (cid:48) ν ( q ).” In all three panels, wesee that the “isotropic S (cid:48) ν ( q )” curves are far from the fullresults (solid red curves), meaning that the anisotropy in S (cid:48) ν ( q ) plays an important role in determining the totalmodulation pattern. We find the same conclusion for theother materials and for other m χ , ω min values.We can further dissect the anisotropy in S (cid:48) ν ( q ) by com-puting the daily modulation with ω − ν, k or Y j · (cid:15) ν, k ,j madeisotropic by the same prescription as above; these are la-beled “isotropic ω − ν, k ” (dotted blue curves) and “isotropic Y j · (cid:15) ν, k ,j ” (dashed green curves) in Fig. 2, respectively.We see that the anisotropy in the Y j · (cid:15) ν, k ,j factor con-tributes the most to daily modulation, as making itisotropic leads to the most significant deviations from thefull results. We find the same is true for other materials. - - - Figure 3.
Left:
Daily modulation for a h-BN target with various experimental thresholds, ω min , assuming dark photonmediated scattering and m χ = 100 keV. Right:
Differential rate at t = 0 for the same process assuming σ e = 10 − cm .The daily modulation pattern is drastically different depending on whether the optical phonon modes just below 100 meV areincluded or excluded. [ h ] Figure 4. Effect of the crystal target orientation on the daily modulation pattern, for a sapphire target and the light darkphoton mediator model as an example. The default orientation is the one adopted in Refs. [48, 49, 63] for which v e ( t ) is givenby Eq. (10), and the alternative orientation is achieved by rotating the crystal z axis by 60 ◦ clockwise around ˆ n = ( ˆ x + ˆ y + ˆ z ) / √ We have also examined the effect of setting e i G · x j → S (cid:48) ν ( q ) while leaving both ω − ν, k and Y j · (cid:15) ν, k ,j intact.This has a visible impact only when the region of q space just outside the 1BZ has a significant contribu-tion to the rate; as q moves farther away from the 1BZ,summing over contributions from many different G vec-tors mitigates the effect. For the dark photon mediatormodel, this explains the enhanced daily modulation at m χ (cid:38) MeV (see Fig. 5 below). For the light hadrophilicscalar mediator model, there is no significant effect sincethe q integral is dominated by small q . For the heavyhadrophilic scalar mediator model, in contrast, the q in-tegral is dominated by large q , so the enhancement hap-pens in a window around m χ ∼ MeV (see Fig. 7 below).
3. Effects of Experimental Setup
The daily modulation pattern can also be significantlyaffected by experimental factors, including in particularthe detector’s energy threshold and the orientation of thetarget crystal. The energy threshold ω min can be impor-tant if phonon modes at different energies have differentmodulation patterns. As an example, we show in the leftpanel of Fig. 3 the daily modulation in h-BN for severaldifferent values of ω min , for the dark photon mediatormodel with m χ = 100 keV. The distinct daily modula-tion curves can be understood from the differential rateplot in the right panel of Fig. 3. We see that the phononmodes just below 100 meV dominate the total rate, solong as ω min is below this, and they drive the daily mod- - - - - - - - - - - - - - - - - - - - - - -
10 25 50 75 100 250 500 10 - -
10 25 50 75 100 250 500 10 - - Figure 5.
Top:
Projected reach for the dark photon mediator model assuming 1 meV and 20 meV energy thresholds andone kg-year exposure. Solid curves show the 95% confidence level (CL) exclusion limits in the case of zero observed events,assuming no background. Dashed curves and the associated ± σ bands show the modulation reach for DM masses with morethan 1% daily modulation, i.e. cross sections for which we can reject the non-modulating hypothesis and establish the statisticalsignificance of a modulating signal, as explained in App. B. Bottom:
Daily modulation amplitudes f mod , defined in Eq. (15),for the same energy thresholds for selected DM masses. Results are shown only for m χ values where a material has substantialreach and f mod > − . ulation pattern. On the other hand, if ω min >
100 meV,these modes are no longer accessible, and the daily mod-ulation is instead induced by phonon modes at energieshigher than about 175 meV, for which the rate has a verydifferent time dependence.Meanwhile, the orientation of the crystal determinesthe function v e ( t ), and hence the daily modulation pat-tern. As an example, Fig. 4 compares the daily mod-ulation patterns between our default setup, given inEq. (10), and an (arbitrarily chosen) alternative orienta-tion where the crystal z axis is rotated by 60 ◦ clockwisearound ˆ n = ( ˆ x + ˆ y + ˆ z ) / √ III. TARGET COMPARISON
Having discussed the physics underlying daily modula-tion, we now consider concrete target materials. Amongthe 26 materials studied [66], 19 are observed to havemore than 1% daily modulation for some DM masses in at least one of the benchmark models considered. In thissection, we focus on the following five: Al O , SiO , SiC,CaWO and h-BN. Among them, Al O , SiO and SiChave been proposed and recommended for near-futurephonon-based experiments, Al O and CaWO are inuse in the CRESST experiment, while h-BN is a highlyanisotropic target that we have found to have exception-ally large daily modulation. We supplement this analysiswith the remaining 14 materials with more than 1% dailymodulation (AlN, CaF , GaN, GaSb, InSb, LiF, MgF ,MgO, NaF, PbS, PbSe, PbTe, ZnO, ZnS) in App. A.Our main results are shown in Figs. 5, 6 and 7, forthe dark photon mediator model and the light and heavyhadrophilic scalar mediator models, respectively. In thetop panels of each figure, we show both the projectedexclusion limits (solid) and the cross sections needed todistinguish the modulating signal and a non-modulatinghypothesis in the event of discovery (dashed and shaded ± σ bands), assuming 1 and 20 meV energy thresholds.For the solid curves, we set t = 0 when computing the - - - - - - - - - - - - - - - - - - - - - -
10 25 50 75 100 250 500 10 - -
10 25 50 75 100 250 500 10 - - Figure 6. Same as Fig. 5, for the light hadrophilic scalar mediator model. rates for concreteness, and assume 3 events per kilogram-year exposure, corresponding to 95 % confidence level(CL) exclusion in a background-free experiment. The re-sults for Al O , CaWO and SiO were computed previ-ously in Ref. [50] (numerical errors in some of the materi-als in early versions of that reference have been correctedhere and on the interactive webpage [66]), and here weperform the calculation also for SiC and h-BN. For thedashed curves and the shaded bands for the modulationreach, we compute the number of events needed to rejectthe constant rate hypothesis at the 95 % confidence levelby a prescription discussed in App. B; they are truncatedwhere the daily modulation falls below 1%.In the lower panels of Figs. 5, 6 and 7, we quantifythe amount of daily modulation for several representativeDM masses by f mod ≡ R max − R min (cid:104) R (cid:105) , (15)where R max(min) is the maximum (minimum) direct de-tection rate throughout the day, and (cid:104) R (cid:105) is the dailyaverage. We shall refer to f mod as the daily modulationamplitude. The f mod plots give us an overview of theamount of daily modulation to expect. More detailed in-formation on the daily modulation signal can be gained by plotting R ( t ) / (cid:104) R (cid:105) , as in Figs. 2, 3 and 4, for each DMmass and energy threshold; we provide these plots on theinteractive webpage [66].We have considered detector energy thresholds ω min =1 meV and 20 meV. For the dark photon mediator model(Fig. 5), the energy threshold does not have a signifi-cant impact on either the reach or the daily modulationamplitude, except at the lowest m χ values. This is be-cause gapped optical phonons dominate the rate as longas they are above ω min and the DM is heavy enough toexcite them. For the hadrophilic scalar mediator mod-els (Figs. 6 and 7), on the other hand, gapless acousticphonons dominate and, as a result, both the reach andthe daily modulation amplitude are sensitive to ω min .Generally, a higher energy threshold tends to amplifythe daily modulation since the kinematically accessiblephase space becomes limited, as discussed in detail inSec. II B 1. Similarly, the daily modulation amplitudetends to increase at the lowest m χ considered because ofphase space restrictions. The enhanced daily modulationin these cases comes at the price of a lower total rate, sothere is a trade-off between better overall sensitivity anda higher daily modulation signal. This is reflected by thedashed modulation reach curves in the top panels of eachfigure, which ascend at lower masses since the rate alsovanishes. - - - - - - - - - - - - - - - - - -
10 25 50 75 100 250 500 10 - -
10 25 50 75 100 250 500 10 - - Figure 7. Same as Fig. 5, for the heavy hadrophilic scalar mediator model.
From Figs. 5, 6 and 7, we see that h-BN consistentlyoutperforms all other materials in terms of the daily mod-ulation amplitude, which reaches O (1) for some m χ and ω min values. This is due to the layered crystal structurewhich means that the momentum transfers perpendicu-lar and parallel to the layers lead to very different targetresponses. Among the other materials, Al O , CaWO and SiC are also competitive targets for the dark photonmediator model at m χ (cid:46)
100 keV, and CaWO showspercent level daily modulation across a wide range of DMmasses for the heavy scalar mediator model.It is also worth noting that the modulation reachcurves and f mod often exhibit a nontrivial dependenceon m χ . In particular, for given target material and ω min ,there can be m χ values where the modulation signal di-minishes. For example, for dark photon mediated scat-tering, h-BN with ω min = 1 meV has two such low- f mod mass points at around 20 keV and 200 keV, correspond-ing to the peaks of the modulation reach curve in thetop-left panel of Fig. 5.Generally, low- f mod points at low m χ result from thechange in ˆ q · ˆ v e favored by the kinematic function. As dis-cussed in Sec. II B 1, as m χ increases, the favored ˆ q · ˆ v e increases from − ˆ v e changes with time( e.g. as in Eq. (10)), a given ˆ q · ˆ v e probes the crystal’s S (cid:48) ν ( q ) along a set of ˆ q directions that modulates, and the modulation pattern depends on the kinematically fa-vored ˆ q · ˆ v e value. We verify this expectation in the leftpanel of Fig. 8 for h-BN. In this case, the modulationpattern flips as m χ increases from 10 keV to 40 keV, andan approximate cancelation occurs around 20 keV. Note,however, that the daily modulation sensitivity may berecovered by analyzing the differential rates dR ( t ) dω .The low- f mod points at higher m χ , on the other hand,are explained by new phonon modes with different modu-lation patterns becoming kinematically accessible as m χ increases. Again focusing on h-BN as an example, we seefrom the right panel of Fig. 8 that while the dominantphonon modes are the ∼
100 meV modes for m χ = 50 keVand 100 keV, the modes above 150 meV take over as m χ increases to 250 keV. The reduced modulation sensitivityat m χ (cid:39)
200 keV results from the transition between thetwo regimes.
IV. CONCLUSIONS
As new experiments focused on light DM detec-tion with optical and acoustic phonons begin an R&Dphase [61], it is important and timely to understandwhich target crystals have the optimal sensitivity to well-motivated DM models. This includes not only the sensi-0 - - - Figure 8.
Left:
Daily modulation for an h-BN target with various DM masses, assuming dark photon mediated scattering and ω min = 1 meV. The change in modulation pattern is a result of the kinematically favored ˆ q · ˆ v e increasing from − m χ increases. During the transition between different modulation patterns, an intermediate mass value around 20 keV featuresa reduced modulation amplitude, which explains the peak in the modulation reach curve in the top-left panel of Fig. 5. Asimilar effect is also observed for the hadrophilic scalar mediator models in Figs. 6 and 7. Right:
Differential rates at t = 0for several higher m χ assuming σ e = 10 − cm . Another transition between modulation patterns occurs when new phononmodes become dominant as m χ increases, resulting in a second reduced modulation mass point, around 200 keV, in Fig. 5. tivity to the smallest interaction cross section for a givenDM model, but also the ability to extract a smoking gunsignature for DM that can be distinguished from back-ground. Daily modulation provides such a unique finger-print. In this work, we have carried out a comparativestudy of daily modulation signals for several benchmarkmodels, where DM scattering is mediated by a dark pho-ton or hadrophilic scalar mediator. Our results supple-ment the information on the cross section reach obtainedpreviously in Ref. [50], and provide further theoreticalguidance to the optimization of near future phonon-basedexperiments.Our analysis shows that there is often a trade-off be-tween detection rate, modulation amplitude, and exper-imental feasibility. For example, for dark photon me-diated scattering, Al O (sapphire), CaWO and SiO ( α -quartz) outperform h-BN in terms of their sensitiv-ities to the total rate; h-BN’s daily modulation signal,however, is significantly stronger. Still, despite havingthe largest daily modulation amplitude, h-BN will likelybe difficult to fabricate as a large ultra-pure single crystaltarget. Overall, Al O and CaWO provide perhaps theoptimal balance between the overall reach and the dailymodulation signal, and have both already been used indirect detection experiments.Beyond the results presented in this paper, we alsopublish an interactive webpage [66], where additional re-sults can be generated from our calculations of singlephonon excitation rates and their daily modulation. Acknowledgments.
We thank Sin´ead Griffin and KatieInzani for previous collaboration on DFT calculations(published in Ref. [50]) utilized in this work. We alsothank Matt Pyle for useful discussions. This material is based upon work supported by the U.S. Department ofEnergy, Office of Science, Office of High Energy Physics,under Award Number DE-SC0021431, and the QuantumInformation Science Enabled Discovery (QuantISED) forHigh Energy Physics (KA2401032).
Appendix A: Daily Modulation Amplitudes forAdditional Materials
In addition to the materials discussed in the main text,we have also investigated the daily modulation of the fulllist of materials considered in Ref. [50]. Their projectedreach curves were already computed in Ref. [50] and areincluded on the interactive webpage [66]. Here we onlyshow the daily modulation amplitude f mod , defined inEq. (15), as in the lower panels of Figs. 5, 6 and 7 inthe main text. The results for the dark photon mediatormodel, the light hadrophilic scalar mediator model andthe heavy hadrophilic scalar mediator model are shownin Figs. 9, 10 and 11, respectively. Only materials with f mod ≥ − for at least one m χ value are shown in eachcase. Appendix B: Calculation of the Modulation Reach
To establish the statistical significance of a modulatingsignal, we find the expected number of events needed toreject the non-modulating hypothesis using the followingprocedure. For a given DM model, with the DM mass m χ and experimental energy threshold ω min specified, we1
10 25 50 75 100 250 500 10 - - Figure 9. Daily modulation amplitudes for the dark photon mediator model for selected DM masses. Solid and dashed linesassume energy thresholds of 1 meV and 20 meV, respectively. Among the materials studied, only those that have a modulationamplitude greater than 1% for at least one m χ value are shown. As in the lower panels of Figs. 5, 6 and 7 in the main text,the low mass values where the rate diminishes are excluded for each material. Therefore the shown modulation amplitudescorrespond to the mass values where the materials have reach.
10 25 50 75 100 250 500 10 - -
10 25 50 75 100 250 500 10 - - Figure 10. Same as Fig. 9, for the light hadrophilic scalar mediator model. first obtain the modulating signal shape r ( t ) ≡ R ( t ) / (cid:104) R (cid:105) as explained in the main text. We divide a day into N bins = 24 equal-size bins, and denote the bin bound-aries by t k = ( k/N bins ) days. Given an expected numberof events N exp , we simulate a DM signal sample and anon-modulating sample by generating events following aPoisson distribution in each bin, with mean (cid:104) N k (cid:105) sig ≡ N exp (cid:82) t k t k − r ( t ) dt/ day and (cid:104) N k (cid:105) non-mod ≡ N exp /N bins forthe k th bin, respectively. We define our test statisticto be the difference between the Pearson’s χ valueswhen fitting the simulated data to the non-modulatingvs. modulating signal shapes. Concretely, suppose thenumber of events in the k th bin is N k . The test statistic is given byTS = (cid:88) k ( N k − (cid:104) N k (cid:105) non-mod ) (cid:104) N k (cid:105) non-mod − (cid:88) k ( N k − (cid:104) N k (cid:105) sig ) (cid:104) N k (cid:105) sig . (B1)Given N exp , we simulate events according to the mod-ulating (DM signal) and non-modulating hypotheses for N sample = 10 times each, and obtain the distributionof TS for the modulating and non-modulating samples.For the non-modulating sample, we compute the 95 per-centile value TS non-mod, 95% . For the modulating sig-nal sample, we compute the mean TS sig, mean and the(50 ±
34) percentiles TS sig, ± σ . These numbers tell us towhat extent we can reject the non-modulating hypoth-2
10 25 50 75 100 250 500 10 - -
10 25 50 75 100 250 500 10 - - Figure 11. Same as Fig. 9, for the heavy hadrophilic scalar mediator model. esis: TS non-mod, 95% < TS sig, mean means we can rejectthe non-modulating hypothesis at 95% CL on average,while TS non-mod, 95% < TS sig, ± σ means we can rejectthe non-modulating hypothesis at 95% CL given a ± σ statistical fluctuation of the signal. Repeating the calcu-lation for many values of N exp , we obtain the interpolat-ing functions TS non-mod, 95% ( N exp ), TS sig, mean ( N exp ) andTS sig, ± σ ( N exp ). These allow us to solve for the N exp needed for TS non-mod, 95% to drop below TS sig, mean , and for it to go below TS sig, ± σ . These then translate intocross sections assuming 1 kg-yr exposure, represented bythe modulation reach curves in Figs. 5, 6 and 7. Notethat the procedure here largely follows that in Ref. [48],but we have adopted a different test statistic that we findsimpler to compute and interpret. We have checked thatusing instead the test statistic in Ref. [48] produces verysimilar results in most cases. [1] J. Amar´e et al. , J. Phys. Conf. Ser. , 012014 (2020),arXiv:1910.13365 [astro-ph.IM].[2] F. Pr¨obst et al. , Proceedings, 7th International Work-shop on Topics in Astroparticle and Underground Physics(TAUP 2001): Gran Sasso, Assergi, L’Aquila, Italy,September 8-12, 2001 , Nucl. Phys. Proc. Suppl. , 67(2002).[3] G. Angloher et al. (CRESST), Eur. Phys. J. C , 43(2019), arXiv:1809.03753 [hep-ph].[4] H. Kluck et al. (CRESST), J. Phys. Conf. Ser. ,012038 (2020).[5] R. Bernabei et al. , Prog. Part. Nucl. Phys. , 103810(2020).[6] A. Aguilar-Arevalo et al. (DAMIC), Phys. Rev. Lett. , 241803 (2020), arXiv:2007.15622 [astro-ph.CO].[7] P. Agnes et al. (DarkSide), Phys. Rev. Lett. , 081307(2018), arXiv:1802.06994 [astro-ph.HE].[8] J. H. Jo (DM-Ice), Proceedings, 38th International Con-ference on High Energy Physics (ICHEP 2016): Chicago,IL, USA, August 3-10, 2016 , PoS
ICHEP2016 , 1223(2017), arXiv:1612.07426 [physics.ins-det].[9] K. Kim (KIMS), in
Proceedings, Meeting of the APS Di-vision of Particles and Fields (DPF 2015): Ann Arbor,Michigan, USA, 4-8 Aug 2015 (2015) arXiv:1511.00023[physics.ins-det].[10] D. S. Akerib et al. (LUX), Phys. Rev. D , 042001(2020), arXiv:1907.06272 [astro-ph.CO]. [11] D. S. Akerib et al. (LUX), (2020), arXiv:2003.11141[astro-ph.CO].[12] M. Antonello et al. (SABRE), Eur. Phys. J. C , 363(2019), arXiv:1806.09340 [physics.ins-det].[13] R. Agnese et al. (SuperCDMS), Phys. Rev. D , 062001(2019), arXiv:1808.09098 [astro-ph.CO].[14] I. Alkhatib et al. (SuperCDMS), (2020),arXiv:2007.14289 [hep-ex].[15] E. Aprile et al. (XENON), Phys. Rev. Lett. , 251801(2019), arXiv:1907.11485 [hep-ex].[16] L. J. Hall, K. Jedamzik, J. March-Russell, and S. M.West, JHEP , 080 (2010), arXiv:0911.1120 [hep-ph].[17] N. Bernal, M. Heikinheimo, T. Tenkanen, K. Tuomi-nen, and V. Vaskonen, Int. J. Mod. Phys. A , 1730023(2017), arXiv:1706.07442 [hep-ph].[18] M. J. Strassler and K. M. Zurek, Phys. Lett. B , 374(2007), arXiv:hep-ph/0604261.[19] N. Arkani-Hamed and N. Weiner, JHEP , 104 (2008),arXiv:0810.0714 [hep-ph].[20] C. Cheung, J. T. Ruderman, L.-T. Wang, and I. Yavin,Phys. Rev. D80 , 035008 (2009), arXiv:0902.3246 [hep-ph].[21] D. E. Morrissey, D. Poland, and K. M. Zurek, JHEP ,050 (2009), arXiv:0904.2567 [hep-ph].[22] D. E. Kaplan, M. A. Luty, and K. M. Zurek, Phys. Rev. D79 , 115016 (2009), arXiv:0901.4117 [hep-ph]. [23] T. Cohen, D. J. Phalen, A. Pierce, and K. M. Zurek,Phys. Rev. D82 , 056001 (2010), arXiv:1005.1655 [hep-ph].[24] K. Petraki and R. R. Volkas, Int. J. Mod. Phys. A ,1330028 (2013), arXiv:1305.4939 [hep-ph].[25] K. M. Zurek, Phys. Rept. , 91 (2014),arXiv:1308.0338 [hep-ph].[26] Y. Hochberg, E. Kuflik, T. Volansky, and J. G. Wacker,Phys. Rev. Lett. , 171301 (2014), arXiv:1402.5143[hep-ph].[27] Y. Hochberg, E. Kuflik, H. Murayama, T. Volansky,and J. G. Wacker, Phys. Rev. Lett. , 021301 (2015),arXiv:1411.3727 [hep-ph].[28] R. Essig, J. Mardon, and T. Volansky, Phys. Rev. D85 ,076007 (2012), arXiv:1108.5383 [hep-ph].[29] P. W. Graham, D. E. Kaplan, S. Rajendran, and M. T.Walters, Phys. Dark Univ. , 32 (2012), arXiv:1203.2531[hep-ph].[30] S. K. Lee, M. Lisanti, S. Mishra-Sharma, and B. R.Safdi, Phys. Rev. D92 , 083517 (2015), arXiv:1508.07361[hep-ph].[31] R. Essig, M. Fernandez-Serra, J. Mardon, A. Soto,T. Volansky, and T.-T. Yu, JHEP , 046 (2016),arXiv:1509.01598 [hep-ph].[32] S. Derenzo, R. Essig, A. Massari, A. Soto, and T.-T.Yu, Phys. Rev. D96 , 016026 (2017), arXiv:1607.01009[hep-ph].[33] Y. Hochberg, Y. Kahn, M. Lisanti, C. G. Tully,and K. M. Zurek, Phys. Lett.
B772 , 239 (2017),arXiv:1606.08849 [hep-ph].[34] R. Essig, T. Volansky, and T.-T. Yu, Phys. Rev.
D96 ,043017 (2017), arXiv:1703.00910 [hep-ph].[35] N. A. Kurinsky, T. C. Yu, Y. Hochberg, and B. Cabrera,Phys. Rev.
D99 , 123005 (2019), arXiv:1901.07569 [hep-ex].[36] C. Blanco, J. I. Collar, Y. Kahn, and B. Lillard, Phys.Rev. D , 056001 (2020), arXiv:1912.02822 [hep-ph].[37] R. Catena, T. Emken, N. A. Spaldin, and W. Tarantino,Phys. Rev. Res. , 033195 (2020), arXiv:1912.08204 [hep-ph].[38] M. Settimo (DAMIC, DAMIC-M), in (2020) arXiv:2003.09497 [hep-ex].[39] L. Barak et al. (SENSEI), Phys. Rev. Lett. , 171802(2020), arXiv:2004.11378 [astro-ph.CO].[40] D. W. Amaral et al. (SuperCDMS), Phys. Rev. D ,091101 (2020), arXiv:2005.14067 [hep-ex].[41] K. Schutz and K. M. Zurek, Phys. Rev. Lett. , 121302(2016), arXiv:1604.08206 [hep-ph].[42] S. Knapen, T. Lin, and K. M. Zurek, Phys. Rev. D95 ,056019 (2017), arXiv:1611.06228 [hep-ph]. [43] F. Acanfora, A. Esposito, and A. D. Polosa, Eur. Phys.J.
C79 , 549 (2019), arXiv:1902.02361 [hep-ph].[44] A. Caputo, A. Esposito, and A. D. Polosa, Phys. Rev.D , 116007 (2019), arXiv:1907.10635 [hep-ph].[45] G. Baym, D. Beck, J. P. Filippini, C. Pethick,and J. Shelton, Phys. Rev. D , 035014 (2020),arXiv:2005.08824 [hep-ph].[46] A. Caputo, A. Esposito, F. Piccinini, A. D. Polosa, andG. Rossi, (2020), arXiv:2012.01432 [hep-ph].[47] S. Knapen, T. Lin, M. Pyle, and K. M. Zurek, Phys.Lett.
B785 , 386 (2018), arXiv:1712.06598 [hep-ph].[48] S. Griffin, S. Knapen, T. Lin, and K. M. Zurek, Phys.Rev.
D98 , 115034 (2018), arXiv:1807.10291 [hep-ph].[49] T. Trickle, Z. Zhang, K. M. Zurek, K. Inzani, and S. Grif-fin, JHEP , 036 (2020), arXiv:1910.08092 [hep-ph].[50] S. M. Griffin, K. Inzani, T. Trickle, Z. Zhang, andK. M. Zurek, Phys. Rev. D , 055004 (2020),arXiv:1910.10716 [hep-ph].[51] B. Campbell-Deem, P. Cox, S. Knapen, T. Lin, andT. Melia, Phys. Rev. D , 036006 (2020), [Erratum:Phys.Rev.D 102, 019904 (2020)], arXiv:1911.03482 [hep-ph].[52] S. M. Griffin, Y. Hochberg, K. Inzani, N. Kurinsky,T. Lin, and T. C. Yu, (2020), arXiv:2008.08560 [hep-ph].[53] T. Trickle, Z. Zhang, and K. M. Zurek, (2020),arXiv:2009.13534 [hep-ph].[54] Y. Kahn, G. Krnjaic, and B. Mandava, (2020),arXiv:2011.09477 [hep-ph].[55] S. Knapen, J. Kozaczuk, and T. Lin, (2020),arXiv:2011.09496 [hep-ph].[56] M. Pyle, E. Feliciano-Figueroa, and B. Sadoulet, (2015),arXiv:1503.01200 [astro-ph.IM].[57] H. J. Maris, G. M. Seidel, and D. Stein, Phys. Rev. Lett. , 181303 (2017), arXiv:1706.00117 [astro-ph.IM].[58] J. Rothe et al. , J. Low Temp. Phys. , 1160 (2018).[59] I. Colantoni et al. , J. Low Temp. Phys. , 593 (2020).[60] C. Fink et al. , AIP Adv. , 085221 (2020),arXiv:2004.10257 [physics.ins-det].[61] C. Chang et al. , “Snowmass 2021 Letter of Interest: TheTESSARACT Dark Matter Project,” (2020).[62] A. K. Drukier, K. Freese, and D. N. Spergel, Phys. Rev.D , 3495 (1986).[63] A. Coskuner, A. Mitridate, A. Olivares, and K. M.Zurek, (2019), arXiv:1909.09170 [hep-ph].[64] R. M. Geilhufe, F. Kahlhoefer, and M. W. Winkler,(2019), arXiv:1910.02091 [hep-ph].[65] https://github.com/tanner-trickle/dm-phonon-scatter (cid:135) .[66] https://demo-phonon-web-app.herokuapp.com .[67] A. Togo, “Phonon Database at Kyoto University,” http://phonondb.mtl.kyoto-u.ac.jp .[68] A. Togo and I. Tanaka, Scripta Materialia , 1 (2015).[69] R. H. Helm, Phys. Rev.104