Dark Matter Daily Modulation With Anisotropic Organic Crystals

Aromatic organic compounds, because of their small excitation energies ~ O(few eV) and scintillating properties, are promising targets for detecting dark matter of mass ~ O(few MeV). Additionally, their planar molecular structures lead to large anisotropies in the electronic wavefunctions, yielding a significant daily modulation in the event rate expected to be observed in crystals of these molecules. We characterize the daily modulation rate of dark matter interacting with an anisotropic scintillating organic crystal such as trans-stilbene, and show that daily modulation is an ~ O(1) fraction of the total rate for small DM masses and comparable to, or larger than, the ~ 10% annual modulation fraction at large DM masses. As we discuss in detail, this modulation provides significant leverage for detecting or excluding dark matter scattering, even in the presence of a non-negligible background rate. Assuming a non-modulating background rate of 1/min/kg that scales with total exposure, we find that a 100 kg yr experiment is sensitive to the cross section corresponding to the correct relic density for dark matter masses between 1.3-14 MeV (1.5-1000 MeV) if dark matter interacts via a heavy (light) mediator. This modulation can be understood using an effective velocity scale v* = Delta E/q*, where Delta E is the electronic transition energy and q* is a characteristic momentum scale of the electronic orbitals. We also characterize promising future directions for development of scintillating organic crystals as dark matter detectors.


I. INTRODUCTION
Dark matter-electron scattering is a promising search strategy for sub-GeV dark matter (DM) . In molecular or solid-state systems, atoms are close enough that electronic orbitals overlap significantly, lowering the electronic excitation energies to the eV scale and thus allowing detection of DM particles with ∼ MeV-scale mass which carry eV-scale kinetic energy. Moreover, solidstate systems can exhibit anisotropic electronic wavefunctions (see for example [9,16,20,27,28,33]), enabling directional detection schemes which leverage the characteristic signature of the daily modulation of the direction of the DM wind in the lab frame (first noted in the context of multiple scattering from terrestrial overburden in [54][55][56], followed by the connection to directional detection in [57]).
In this paper, we focus on, and advocate for, a particular class of detector materials for DM-electron scattering: aromatic organic crystals. These compounds have numerous advantages, both practical and theoretical, over existing detectors. Their molecular structures consist primarily of hexagonal carbon rings with alternating single and double bonds (see Fig. 1), and the excited molecular electronic levels at O(5 eV) above the ground state can de-excite by emitting a scintillation photon with O(1) bulk quantum efficiency. As a subset of the authors showed in a previous paper [30], simple organic liquids like benzene and its analogs have superior reach per unit target mass compared to single-electron threshold semiconductor detectors for DM heavier than about 10 MeV. Furthermore, as noted in [30] and as discussed at length in the present work, the planar structure of the molecules leads to a marked anisotropy in the electronic wavefunctions which is absent in silicon and germanium detectors and which allows for directional detection: even though the scintillation emission is isotropic, the excitation rate to the scintillation level depends strongly on the direction of the incoming DM.
Building on our previous work, in this paper we focus on larger organic molecules which are solids at room temperature with known bulk quantum efficiencies at cryogenic temperatures. Furthermore, we focus on molecules with reduced in-plane symmetry. This additional asymmetry means that the lowest-energy transitions are not suppressed, distinguishing them from simpler molecules like benzene. The crystal lattice effects in these organic crystals are small enough that the electronic structure closely resembles that of the isolated molecules. Singlecrystal scintillators can have Avogadro's number of unit cells containing the same relative orientations of the molecules, allowing the anisotropic response to persist. From a practical standpoint, single-crystal samples of trans-stilbene (t-stilbene), which we focus on in this paper, can be manufactured at kilogram scale with order-1 scintillation efficiency, such that even a single scintillation photon produced from a DM scattering event has a high probability of being detected from a large-mass 2 sample [58].
These practicalities suggest that running a t-stilbene experiment with a large exposure is feasible in the near future, so we consider possible interpretations of plausible near-future experimental data. As most sources of background noise (such as radioactive impurities in the target material) are constant in time, daily modulation provides a way to detect dark matter even without reducing the background rate to zero. Cosmic rays, the primary time-varying external background, have a daily modulation that has been constrained to be below the level of 0.1% in underground facilities [59], so a modulating signal with a significantly larger amplitude would provide a clear detection of dark matter.
For example, assuming an average observed rate of 1/60 Hz kg −1 , and without incorporating expectations for daily modulation, the assumed (constant) background rate would limit the reach to approximately an order of magnitude above the interesting parameter space, with no prospects for improvement over time or with a larger experiment. However, a 1 kg·yr exposure reaches DM relic density targets if the data are interpreted with the expected modulation information. The sensitivity of a modulation analysis continues to scale with (exposure) 1/2 even without mitigating backgrounds, so a larger 100 kg·yr exposure improves the reach by an order of magnitude. This probes the relic density target for DM masses 1.3 m χ 14 MeV if the DM interacts through a heavy mediator, or the range 1.5 m χ 1000 MeV if the DM achieves its relic abundance through freezein via a light kinetically-mixed dark photon. Solid-state organic scintillator detectors would therefore greatly reduce the necessity for a low-threshold zero-background experiment in order to conclusively discover or exclude DM.
As a consequence of our analysis of organic crystals, we point out a generic feature of DM-electron scattering in condensed matter systems: daily modulation is governed by the relationship between the DM velocity and an effective electron velocity v * ≡ ∆E/q * , where ∆E is the energy of an electronic transition and q * is the typical momentum scale for the electron wavefunctions governing the transition. A necessary condition for daily modulation is anisotropy of the molecular form factor for electronic transitions, but for the anisotropy to be kinematically accessible, the transition needs to be either near a kinematic threshold (for small DM masses) or have v * close to, but smaller than, the maximum DM velocity in the lab frame v max (for large DM masses). Intriguingly, this suggests that daily modulation, like annual modulation, is driven by the high-velocity tail of the DM velocity distribution. The centrality of these kinematic relations has recently been noted by [38], which discusses v * in the context of maximizing the total rate, and by [60], which uses a related v * to study modulation in the context of single-phonon production.
This paper is organized as follows. In Sec. II, we review the quantum chemistry relevant for describing the  Fig. 5 it can be seen that the packing of the molecules is not the same for the two columns. The bond lengths and angles in the two independent centrosymmetric molecules are listed in Fig. 6. The standard deviations in the bond lengths and angles as obtained from the variance-covariance matrix of the coordinates (Darlow, 1960) are 0.0015 A_ for C-C and 0.11 ° for C-C-C. Owing to the disorder described above we have multiplied these values by a factor of two for the ~ molecule. Corresponding bond lengths and C-C-C valence angles in the two molecules do not show significant differences. The disorder of the molecules is so small that it does not affect the thermal parameters to a large extent. The value of (U~(principal axis)) x/2 for the C atoms is 0.0217 A 2 for the and 0.0177 A_ 2 for the fl molecules.
Comparison with the earlier structure determinations shows that the present work has a considerably higher accuracy. The two centrosymmetric molecules are approximately planar. The torsion angle ~ around the C-C(phenyl) bond is 3.4 ° for the a and 6.9 ° for the fl molecule. The C6H5-C groups show only small devia-FIG. 1. Top: The chemical structure of trans-stilbene. Following the convention common in organic chemistry, vertices are taken to be carbon atoms, single lines are carboncarbon single bonds, and double lines are carbon-carbon double bonds. Our numbering convention for the atoms is shown at each vertex, along with theL andM unit vectors used in our coordinate system. Bottom: A diagram of the unit cell of the trans-stilbene crystal, adopted from Ref. [62]. Here b is the axis of symmetry for the crystal: the positions of the molecules in the middle row are related to those of the upper or lower rows by a translation of 1 2 b and a rotation of 180 • about the b axis. The long axis of each molecule (Li) is not perfectly perpendicular to the b crystal axis. Molecules (1) and (−1) are shown with thin lines, (2) and (−2) with thick lines. molecular orbitals in t-stilbene. In Sec. III we describe the crystal structure of t-stilbene and justify our use of the isolated-molecule orbitals based on experimental measurements of the absorption and emission spectra. In Sec. IV we set up our calculation of the DM scattering rate, including the relevant molecular form factors and the daily modulation of the velocity distribution. In Sec. V we determine the daily modulation signal from DM in the Standard Halo Model (SHM) and derive a convenient test statistic for daily modulation in the presence of a non-modulating background. In Sec. VI we study the kinematic features of t-stilbene which lead to a large modulation amplitude, and explore how a system with a different v * could lead to large daily modulation even for DM up to the GeV scale. We conclude in Sec. VII.

II. MOLECULAR ORBITAL MODEL
Here we determine the many-electron wavefunctions for the electronic transitions which are responsible for the electronic-optical properties of trans-stilbene, following methods appropriate for all aromatic molecules. As shown in Fig. 1 (top), t-stilbene is an alternant hydrocarbon hosting two phenyl groups joined by an ethene bridge in the trans configuration (trans-1,2-diphenylethylene). This 14-carbon molecule is planar in the solid state and presents C 2h molecular symmetry [63,64], which is a Z 2 × Z 2 symmetry group composed of a two-fold symmetry axis, a center of inversion, a horizontal mirror plane, and the identity. Since this is the symmetry of the Hamiltonian, the electronic states of t-stilbene should transform as irreducible representations of C 2h . In order to construct electronic states which accurately describe the energy eigenstates of t-stilbene, we first construct the Hückel molecular orbitals (HMOs) using a simple Linear Combination of Atomic Orbitals (LCAO) model taking into account only direct bonding interactions. We then take into account the configurational interactions and construct fully antisymmetric many-body states following the method of Pople, Pariser, and Parr (PPP) [65][66][67]. The HMOs, Ψ i , of t-stilbene are constructed as linear combinations of Slater-type atomic orbitals (SAOs) where c j i are the coefficients to be determined, φ 2pz are the atomic orbitals, and R i are the equilibrium locations of the carbon nuclei using the numbering conventions in Fig. 1. The 2p z Slater atomic orbital is parameterized as where a 0 = (αm e ) −1 is the Bohr radius and Z eff = 3.15 is the effective nuclear charge of the carbon 2p z orbital [66]. The HMOs diagonalize the 14 × 14 core Hamiltonian matrix, H lm = φ l |H core |φ m , where r|φ m = φ 2pz (r − R m ) and only bonding atoms interact. Diagonalization is done by solving the following system, The core Hamiltonian contains two types of matrix elements; diagonal on-site energies and off-diagonal interaction energies, with values given in Appendix A. The onsite energy is an empirical quantity which is determined by the atomic species, while the off-diagonal energy is a measure of the bonded nuclear interactions determined by the effective charge and bond length. The ambiguity in the construction of degenerate states is resolved by requiring that all electronic states transform as irreducible representations of the C 2h point group. Since the HMOs are constructed from 14 SAOs, diagonalization of the core Hamiltonian results in 14 HMOs Ψ 1 , . . . , Ψ 14 , numbered in order of increasing energy (up to degeneracies). See [30] for an example of this procedure performed on the simpler benzene molecule containing only six carbon atoms.
In order to form antisymmetric many-electron wave functions, we take Slater determinants of the filled HMOs. The ground state |g is approximately given by the following combination of orbitals, where | · · · | denotes the antisymmetrized product of the HMOs and Ψ i is the opposite spin state as Ψ i . Following standard conventions in quantum chemistry, we label some Z 2 symmetries of the many-electron wavefunction by A(B) and g(u), corresponding to (anti)symmetry with respect to transformation under 180 • rotation about the z-axis normal to the molecular plane and inversion through the center of mass, respectively. Notice that the ground state represents an electronic configuration in which the lowest 7 HMO's are filled by pairs of electrons in the spin-singlet configuration. Therefore, the ground state transforms as the A g representation, being totally symmetric under the transformations in the C 2h group. This is a generic feature of the ground state of alternant hydrocarbons. Reduction of the 14-dimensional t-stilbene representation of C 2h predicts 7 A g and 7 B u states [63]. These states correspond to the multi-electron configurations of the HMOs which each have either A u or B g symmetry [64].
We construct the multi-electron states starting with the A g ground state and proceeding upwards in energy through the one-electron singlet excitation configurations, ψ j i as follows, The electronic repulsion is taken into account by appending the electronic two-body interaction to the core Hückel Hamiltonian: where α is the fine-structure constant, r ij = |r i − r j |, and the sum runs over all pairs of electrons ij with i = j. This PPP Hamiltonian perturbs the energy levels of the ψ j i configurations and mixes degenerate states of like symmetries to produce PPP HMOs of either A g or B u symmetry [61,63,64]. The PPP energy eigenstates are expressed as a linear combination of ψ j i , The leading coefficients d (n) ij (known as configuration amplitudes) for the first n = 1 . . . 8 excited states of t-stilbene are tabulated in Tab. I [61]. To better match our analysis to experimental data, we take the experimentally-determined energy eigenvalues, which are also listed in Tab. I [61].
The spin-singlet configurations we have focused on are responsible for the radiative de-excitations known as fluorescence that could be seen with single-photon detectors in a DM experiment. The triplet states are classically forbidden from decaying down to the ground state and hence are responsible for the delayed fluorescence component of photoluminescence which generically has a significantly lower quantum yield [68].

III. CRYSTAL STRUCTURE
In the solid state, t-stilbene is a monoclinic crystal belonging to the space group C5 2h (P 2 1 /c), with unit cell parameters a = 12.29Å, b = 5.66Å, c = 15.48Å, and γ = 112 • [62]. In this convention, the crystal coordinate basis is defined by a unit vectorb = b/b that is orthogonal to bothâ andĉ, and γ is the opening angle betweenâ andĉ. The coordinate system is right-handed, such that c × a b. To form an orthonormal basis we define an a unit vector,â =b ×ĉ. The crystal is symmetric with respect to translations of a and c, and to the twofold screw action composed of the translation b/2 with 180 • rotation aboutb.
Four distinct molecules of t-stilbene inhabit each unit cell of the crystal (see Fig. 1, bottom): an M 1 and M 2 with different orientations, and an M −1 and M −2 , which are the images of molecules M 1,2 (respectively) under the twofold screw action along b [69]. In Tab. II we provide the orientations of each of the four molecules in terms of their unit vectorsL andM identified in Fig. 1, in a crystal coordinate system where theẑ direction is assigned to the b symmetry axis. The position of each molecule within the crystal is listed in Refs. [62,69], but the DMstilbene scattering rate depends only on their rotational orientation because the kinematics of the scattering process do not permit coherent scattering over an entire unit cell (see Secs. IV and V below). The molecular orbital model derived in the preceding section is empirically valid for macroscopic singlecrystal samples of t-stilbene. Although lattice effects are known to perturb the energy bands of molecular crystals [70], the Davydov splitting of the molecular bands, which is due to the dipole-dipole (and to a lesser extent quadrupole-quadrupole) interaction of neighboring molecules in a molecular crystal, is known to be very small in single-crystal t-stilbene [71]. Furthermore, the lowest two UV-absorption bands (labeled A and B in the literature) of t-stilbene in the liquid state remain relatively untouched when observed in the solid state [72]. The A band is thought to arise from the g → s 1 through g → s 4 transitions while the B band is thought to arise from the g → s 5 transition [63]. These bands are a direct measurement of transitions between the many-electron configurations described in the previous section which represent the non-interacting single-molecule electronic state, usually taken to be most like the molecular environment of the low-temperature liquid or gas state. Since these bands remain the same in the solid states, we conclude that the PPP model accurately describes the molecules of t-stilbene in the solid state where lattice effects are only very weak.
A mole of t-stilbene has a mass 180.24 g and occupies a volume of 185.69 cm 3 . Thus, a kilogram of detector material can be fabricated from a cube of t-stilbene of 10.1 cm per side, or a 1 cm thin sheet of approximately one square foot. These sizes will be convenient to instrument with conventional photodetectors, such as photomultiplier tubes or charge-coupled devices (CCDs).

IV. RATE CALCULATION
The interaction of DM with the electronic system of a molecule produces a detectable signal through a process analogous to photoluminescence. This process can be separated into two stages: electronic excitation from DM-electron scattering followed by radiative deexcitation. Specifically, DM-electron scattering will produce detectable scintillation photons from t-stilbene at a rate given by [30] Here, N A 6.022 × 10 23 is Avogadro's number; m T is the total detector mass; the molar mass of t-stilbene is m (t-stil) mol = 180.25 g; ρ χ and m χ are the DM mass density and mass, respectively; µ χe is the DM-electron re- |M(q 0 )| 2 is a fiducial DMelectron cross section proportional to the free-particle spin-averaged squared matrix element |M| 2 for χ − e scattering, evaluated at q 0 = αm e ; f χ (u) is the DM velocity distribution in the lab frame; F DM (q) is a form factor parameterizing the fundamental DM-electron interactions, which has the limits F DM (q) = 1 for a contact interaction and F DM (q) = (αm e /q) 2 for a long-range interaction; and ∆E(s i ) is the excitation energy for each s i above the ground state as given in Tab. I.
The detector-dependent quantities are the molecular form factor |f g→si (q)| 2 , representing the transition amplitude from the ground state to a singlet excited state s i , and the bulk fluorescence quantum efficiency Φ F B , representing the probability that a molecular excitation will produce a photon through radiative deexcitation that will exit the detector without being absorbed. As is the case with photoluminescence, the emission lines are broadened by vibrational energy sublevels, thermal motion, and lattice effects. The emission spectrum of the single crystal is a continuum of peaks which closely resembles the molecular and micro-crystalline emission spectra but is modified by self-absorption and lattice effects [73]. Here we focus on computing the form factor and the quantum efficiency, which determine the signal rate, leaving a detailed investigation of the emission spectrum (which determines the signal photon wavelength and hence the detection mechanism) for future work. We compute the molecular form factors using the first 8 singlet transitions since these are the transitions responsible for the first three lowest-lying absorption bands of trans-stilbene [63]. Since the probability of interaction is suppressed for higher energy thresholds (see Fig. 11), it is expected that the rate will be driven primarily by the strongest low-lying excitations. Using a simple, spherically-symmetric DM velocity distribution, the SHM, we also calculate daily modulation effects due to the rotation of the Earth.

A. Molecular Form Factors
The molecular form factor as calculated for PPP configurational states s n is given by the following, where r m is the position of electron m. In the second line we have isolated the contribution from the singlet states which only contain a single-electron excitation above the ground state, and in the third line we have transformed to the basis of HMOs, where the factor of √ 2 is effectively a spin degeneracy factor. The matrix element of HMOs may be readily computed from Eq. (1) in momentum space, where the wavefunctions are simply the momentum-space 2p z orbitals times a product of phase factors determined by the positions of the carbon atoms [30]. An example of the single-molecule form factor for the g → s 1 transition is shown in Fig. 2 (see also App. A), showing strong damping of the form factor beyond a characteristic momentum scale given by q * Z eff /(2a 0 ) 6 keV. There are also secondary "inner" peaks at q 2π/ 1.2 keV where ≈ 0.83 nm is the length scale of the long axis of t-stilbene.
Measurements of the spectrum of trans-stilbene observe three absorption bands, A, B and C, which have been identified primarily with the s 1 , s 5 and s 8 molecular transitions, respectively [63]. Compared with the B and C bands, the A band is larger in magnitude and broader in frequency space, overlapping with the s 2...4 transitions. Our analysis for dark matter scattering reproduces these features: for example, the s 1 transition dominates the scattering rate, comprising 50-70% of the rate both near-threshold and at large m χ . If m χ is large enough that the higher-energy excitations are kinematically accessible, the scattering rate receives secondary contributions from s 3 and s 4 , and typically smaller contributions from s 2 , s 5 and s 8 . The s 6 and s 7 transition rates remain negligibly small at all values of m χ . We provide additional detail regarding the separate molecular transitions that contribute to the scattering rate in Appendix A 2.
In a t-stilbene crystal, the unit cell contains four molecules of different orientations, as described in Section III. As discussed in more detail below, the momentum transfers required to deposit energy above ∆E(s 1 ) are sufficient to localize the interaction to a single molecule within a unit cell, so to compute the rate we treat the scattering as incoherent between different molecules, and we sum over four different squared form factors |f (i) g→s (q)| 2 rotated to give the appropriate orientations of each molecule with respect to q. The interaction rate of DM with a crystal of t-stilbene then scales like the product of the rate calculated as described with the total number of unit cells N uc = N mol /4 in the entire crystal, which is proportional to the total crystal mass m T .

B. Quantum Efficiency
Given the molecular fluorescence quantum efficiency (Φ F ), defined as the ratio of emitted photons to absorbed, the probability Φ F B of a photon exiting the bulk target after an excitation is then given by where a xx is the probability of self absorption. At liquid nitrogen temperatures, Φ F 97% [74][75][76], approaching unity at cryogenic temperatures. Furthermore, using the photoluminescent spectra of t-stilbene in the liquid, microcrystalline, and macroscopic single-crystal state as measured via reflection and transmission, Birks et al. conclude that a xx ≈ 0.35 as the continuous reabsorption and emission of the photon gradually red-shifts the radiation into wavelengths to which the bulk crystal is transparent [73]. Thus, the bulk fluorescence quantum efficiency of t-stilbene is at least Φ F B = 0.63 [58,77], likely approaching 0.65 at cryogenic temperatures, though we use the lower value to be conservative. We propose such a detector to be run around 100K in order to maximize bulk quantum efficiency while maintaining a high enough temperature to run CCD based photo-detectors [41].

C. DM Velocity Distribution
We denote the velocity of a particle in the Milky Way reference frame by v. We adopt the SHM ansatz for the bulk of the Milky Way DM distribution. The DM velocity is then distributed according to where the dispersion σ 0 is related to the velocity of the Local Standard of Rest by σ 0 = v 0 / √ 2 and the normalization is Numerically, the Local Standard of Rest has value v 0 220 km/s and the escape velocity is near v esc 544 km/s [78]; the uncertainties on these values are nonzero, but 10%. We use this velocity distribution in order to facilitate comparison with previous studies, but future exploration of the implications of more realistic velocity distributions will be important for interpreting any future experimental results.
To calculate the rate in Eq. (8), we integrate over the velocity measured in the laboratory u, which is related to the velocity in the Milky Way frame by where v ⊕ (t) is the velocity of the Earth as measured in the Milky Way frame. Following the conventions of [27], we use the energy conservation δ-function to resolve the integral over the velocity distribution: where The momentum transfer q = qq has been related to the velocity by the requirement that ∆E = q · u − q 2 /(2m χ ), and all dependence on the local velocity v ⊕ is contained in v − . The rotation of the Earth over a 24-hour period enters the rate by casting the direction of Earth's velocity (more precisely, the velocity vector of a fixed laboratory location on the Earth's surface) as a function of time [27]: sin θ e sin ϑ sin θ e cos θ e (cos ϑ − 1) cos 2 θ e + sin 2 θ e cos ϑ   , (14) where ϑ(t) = 2π× t 24 h , θ e ≈ 42 • , and we have chosen to align the (x, y) plane of the crystal to be perpendicular to the direction of the DM wind at time t = 0 with the initial orientation of the crystal with respect to rotations about theẑ axis given by β. As is clear from Eq. (13), the only aspects of the Earth velocity vector that we need to know when calculating rates in the context of the SHM are the Earth's speed, for which we adopt |v ⊕ | ≡ v ⊕ = 234 km/s, and the angle between the Earth's north pole and its velocity in the Milky Way frame, θ e . Our formalism is easily extended to other velocity distributions by making the substitution in Eq. (13) the mean velocity of the DM distribution as measured in the Milky Way frame. The vector v is zero by definition for the SHM, but would be nonzero for substructure in the form of a stream.
In terms of g 0 , the total (time-dependent) rate per unit mass is where we emphasize that g 0 (q) is implicitly also a function of v ⊕ (t). Eqs. (12)- (15) indicate that, from the perspective of kinematics alone, the largest modulation will occur when v − (q) modulates around v esc . However, as we will see in Secs. V and VI below, the morphology of the molecular form factors will also play a large role in determining the modulation.

V. DAILY MODULATION REACH
With the molecular form factors in hand, we can compute the total DM-induced excitation rate by summing over the eight lowest transitions s i for a given choice of velocity distribution and DM form factor. As we show in Appendix A, the lowest-energy transition g → s 1 dominates both the daily modulation effect and the total average rate. Near the mass threshold, m χ 10 MeV, the g → s 3 transition contributes at the 20% level, with all other transitions contributing less than 10% of the total rate. Above 10 MeV the s 1 transition remains dominant, accounting for about 50% of the total rate for both DM form factors. With F DM = 1 and m χ 100 MeV the s 3 and s 4 transitions contribute equally, at the 15% level. For the same masses and F DM ∝ 1/q 2 , the s 3 transition provides a larger 20% correction to the total rate, compared to less than 15% from s 4 and less than 10% from each of the other excited states. This behavior is distinct from the case of benzene, where the lowest-energy transition is dipole-forbidden [30]. Fig. 3 shows the modulating rate R(t) over a 24-hour period (one sidereal day) for two different alignment angles of the detector crystal, β = 0 • and β = 90 • , normalized by the average scattering rate, We see that the peak-to-trough modulation amplitude is as large as 60% (10%) for a low-mass (high-mass) DM particle interacting via a heavy mediator, climbing to 70% (25%) for a low-mass (high-mass) DM particle interacting via a light mediator. This is on the same scale, or larger than, the annual modulation amplitude for WIMP-nuclear scattering well above threshold [79][80][81][82][83], as well as for DM-electron scattering in semiconductors at high masses [4,34]. Assessing evidence in favor of a signal will be an important step in making a DM discovery, and the daily modulation is an important handle for improving our statistical power. As we discuss in more detail in Appendix B, the statistical significance that we formally assign to a modulating signal is where k labels the data bins, ν m k is the number of expected events in bin k assuming a modulating signal,   expected rate in the modulating and non-modulating scenario, respectively. The values of ∆L are distributed as a χ 2 distribution of the number of additional degrees of freedom needed to characterize the modulating (as opposed to the non-modulating) signal; in the case of two bins, this would be a χ 2 with two degrees of freedom. Although we focus on the two-bin case in the remainder of this analysis, we emphasize that Eq. (17) is appropriate for any binning of data, including an unbinned analysis. We provide more general explorations of this test statistic in Appendix B.
A particularly simple limit of Eq. (17) is one for which we take two bins per day and describe the modulation simply by a single parameter, the integrated modulation fraction f 2 , defined as the fractional difference in inte-grated rate between the two bins, averaged over a day: For a perfectly sinusoidal signal, f 2 equals the peak-totrough amplitude divided by π. Our choice in Eq. (14) to align the crystalline symmetry axis,b, with the lab frame DM wind at t = 0 ensures that the dominant part of the modulation signal has a 24-hour period, with only small contributions from higher harmonics. In this orientation, the integrated modulation amplitude Eq. (18) is maximized by t 0 ≈ 18 hours, based on the results shown in Fig. 3. This observable is particularly well suited for describing the daily modulation, because it is unaffected by the non-modulating background rate and thus does not require any knowledge of the background.
As explored in detail in Appendix B, this simple binning is amenable to analytic results in the large-N limit  [51], XENON 10 [14], and XENON 1T [50]. The dotted and dashed lines show the 90% CL exclusions that can be set from the total number of events, without considering modulation effects, for R = 1/60 Hz kg −1 (Nevents ≈ 5.26 × 10 5 ) and for Nevents = 0, respectively. The orange shaded regions indicate parameter space that leads to a sufficiently large modulation signal that a 1 kg · year experiment could observe a 3σ detection, given a total observed rate of R = 1/60 Hz kg −1 . The solid black "∆N = 0" lines show the improved limit that can be set from a null result exhibiting no daily modulation but the same total observed rate. Each plot also shows (in blue) a benchmark model from Ref. [13] as a target for the experimental sensitivity. In the FDM = 1 example the scalar DM abundance is set by freeze-out mediated by a dark photon of mass m A = 3mχ, while for FDM ∝ 1/q 2 we show freeze-in via light mediator, m A 3 keV.
of the Skellam distribution or in the small-modulation limit of the Poisson distribution. In each case, we find that the statistical significance we may assign to either the modulating or non-modulating hypothesis based on an experiment in which N tot events are observed is where R is the time average of the signal event rate R(t) defined in Eq. (16), and where T exp is the total exposure time for the experiment. Since the number of signal and background events both grow linearly with exposure, the significance of a modulating signal improves with exposure as long as the integrated modulation fraction f 2 is nonzero. Our Eq. (19) matches the χ 2 sb statistic suggested by Ref. [28].
In Fig. 4 we show the expected results of a 1 kg · year t-stilbene experiment operated under a number of different assumptions. As a benchmark to facilitate comparison with other experiments, we demonstrate the reach with an entirely background-free experiment using no modulation information. The potential for parameter space exclusion in this scenario isσ e 10 −41 cm 2 (few × 10 −41 cm 2 ) for DM interacting with a form factor F DM = 1 (F DM ∝ 1/q 2 ) and with a mass in the range 5 MeV m χ 10 MeV. This is within a factor of 2 or 3 from the expected reach of a silicon CCD experiment like SENSEI or Oscura for an equivalent target mass [3,84]. Taking the more realistic scenario that the observed rate for a 1 kg detector is R = 1 min −1 = 1/60 Hz (including both signal and background components), the future reach depends on analysis strategy. Without leveraging modulation information, the limit we obtain is slightly stronger than the current exclusion from SENSEI [51] at low masses below ∼ 5 MeV, and comparable at higher masses. We also comment in passing on the prospects for the detectability of the (non-modulating) absorption of DM: since ρ T 1 g/cm 3 , we anticipate a rate of ∼ O(1)/kg/min for a dark photon kinetic mixing parameter of 10 −13 , assuming a dielectric loss of order ∼ O(10 −2 ), similar to that in benzene [85] and comparable to those in semiconductors [86]. This setup would set leading limits on dark photons in the mass range 4.2 eV < m A 10 eV. For DM scattering, our sensitivity to exclusion and discovery can be dramatically extended by utilizing the information in the rate via the simple two-bin analysis. Using the significance from Eq. (19), a 1 kg, 1 year t-stilbene experiment that observes a constant R = 1 min −1 = 1/60 Hz event rate can exclude at 90% CL a DM particle with a scattering cross section as small As a demonstration of the utility of daily modulation, we show the 3σ discovery and 90% CL exclusion potential for a trans-stilbene experiment with a background rate of 1/60 Hz kg −1 , for exposures of 0.01 kg · year, 1 kg · year, and 100 kg · year. The dashed lines, labelled "90% CL 1/60 Hz kg −1 ," show the 90% CL exclusion from an analysis that does not consider the daily modulation effects. The inclusion of daily modulation in the statistical analysis allows even the 0.01 kg · year exposure to set a significantly stronger limit onσe. For the FDM = 1 and FDM = (αme/q) 2 form factors we show the benchmark freeze-out and freeze-in models, respectively, from Ref. [13]. asσ e 10 −37 cm 2 . This cross section lies below the well-motivated line from freeze-out production of scalar DM for 2 MeV m χ 7 MeV with a heavy dark photon mediator mediator, and also probes a wide range of masses 2 MeV m χ 200 MeV for freeze-in production through a light mediator [13]. The 3σ discovery reach for a modulating signal for a total R = 1/60 Hz/kg background event rate is nearly as strong, reaching just below (above) the cross sectionσ e 10 −37 cm 2 for F DM = 1 (F DM ∝ 1/q 2 ).
Very meaningfully, as shown in Eq. (19), the discovery or exclusion significance grows with cumulative exposure, even without background mitigation: this improvement in significance is absent in a non-modulating signal. We demonstrate this explicitly in Fig. 5, displaying the 90% CL exclusion and the 3σ discovery reach for a t-stilbene experiment with a constant observed rate R = 1 min −1 kg −1 and increasing exposures of 0.01, 1, 100 kg·yr. For the lowest exposure proposed here, the background rate is very nearly equal to the 2e − rate observed by the SENSEI experiment with a ∼2g detector [51]. The sensitivity improves with √ N tot , so given the assumption of constant total rate in counts per unit time per unit mass, the sensitivity improves with √ exposure. This conservative expectation for scaling of the background rate essentially assumes that bulk events will dominate the background. There will be an irreducible background from the low-energy tail of 14 C decays which would yield only a single scintillation photon, but assuming scintilla-tors can be manufactured with the 10 −18 g/g 14 C levels achieved by Borexino [87,88], the total 14 C decay rate would be 0.01 events/min/kg, well below the background rates we have assumed here. These background rates also include radio contamination from heavy metals (e.g. Th and U) whose beta decay spectrum is not compatible with the one-photon signal that such an experiment will look for. Finally, the cosmic ray background is expected to be the bulk of the exogenous rate and the only major background that might vary over the time scale of a day. This rate can be minimized by running under sufficient overburden, and the daily modulation of this background is constrained to be 0.1% in a deep underground facility [59]. If the background rate were dominated by the dark rate in the photodetector, which would likely scale with area, a large-volume experiment and/or light-focusing scheme would improve the significance even further. We plan to return to these issues in future work.

VI. KINEMATICS AND TARGET SELECTION FOR DAILY MODULATION
Given the large daily modulation amplitude present in t-stilbene, and the associated improvement in discovery and exclusion significance, it is worth examining which characteristics of our target molecule govern the size of the modulation, and whether other choices of organic FIG. 6. Molecular form factors and modulating rates for DM masses near threshold, mχ = 2 MeV. In the contour plots, the gridded shaded regions indicate the kinematically accessible momentum transfers q for the four molecules M1,2,−1,−2 that comprise the unit cell of the crystal, shown at t = 0 and t = 10 h. Here, q is given in the molecular basis, qx = q ·L, qy = q ·M, and the kinematically accessible region is defined by v−(q) < vesc, following Eq. (13). Top left: Contour plot of the molecular form factor |f (q)| 2 for the g → s1 transition in the (qx, qy) plane, with qz = 0. Bottom left: Contour plot for fixed φ ≡ arctan(qy/qx) = 55 • , showing the strong anisotropy in qz with maxima at qz = 0. Top right: The scattering rate (summed over all g → si transitions) as a function of time, R(t), normalized by the average daily rate Ravg. The modulation is dominated by the s1 transition (dashed). Bottom right: A closer look at the form factor near the peak at φ 55 • , plotting |fs 1 | 2 as a function of |q| for fixed θ = 90 • and various φ. molecules could improve the modulation amplitude even further. Indeed, in other systems sensitive to sub-MeV DM (Dirac materials, for example), the modulation amplitude can be even larger, O(1) even for DM masses well above threshold [27,28].

A. Daily modulation in t-stilbene
The peaks of the t-stilbene molecular form factors define a preferred momentum scale q * 6 keV where the rate is largest, so for the purposes of understanding the daily modulation, we may approximate all DM interactions as imparting momentum q * . The s 1 transition has ∆E = 4.2 eV, which defines an effective velocity scale on the same order as v ⊕ 230 km/s. For sufficiently small m χ such that q * /(2m χ ) v esc , Eq. (13) shows that v − (q) will be driven to v esc unless v ⊕ is antiparallel to q, and hence the rate will be nonzero only for a very narrow range of directions of q. Fig. 6 illustrates this phenomenon, with the gridded "bean-shaped" shaded regions representing the kinematically-accessible region v − (q) < v esc overlaid on contour plots of the s 1 molecular form factor at β = 90 • for DM mass of 2 MeV and a heavy-mediator form fac-  Fig. 6 for large DM masses, mχ = 100 MeV. Only the nearly-spherical region near q ∼ 0 with inner boundary qmin 1.6 keV is kinematically forbidden. As a result, the daily modulation amplitude is smaller, driven by the anisotropy of the inner secondary peaks and the tails of the primary peaks.
tor F DM = 1. There are four such regions for the four different t-stilbene orientations within a unit cell. The daily modulation arises from the movement of the kinematically-allowed region over the course of a day, in particular as these regions rotate out of the plane of the molecule and the peaks at q z = 0 become inaccessible.
On the other hand, for sufficiently large m χ , q/(2m χ ) → 0 and v − (q * ) < v esc for any direction ofq. Thus, the kinematically-allowed region in q-space always includes q * but has inner boundary Fig. 7 provides the same information as Fig. 6 except now for a heavier DM particle, with mass m χ = 100 MeV. The form factor remains the same, but the "beans" have now expanded to fill in across the plane, leaving only circular "holes" with inner boundary q min . Because the kinematically-accessible region now includes the full peaks of the form factor, the rate modulation of the course of the day arises only due to the mismatch of the circular inner boundary with the hexagonal symmetry of the form factor and the presence of the inner secondary peaks, compounded by the vector addition ofq andv ⊕ . This leads to a smaller ∼ 10% peak-to-trough modulation amplitude for all m χ 10 MeV.
For DM scattering through a light mediator, F DM = (αm e /q) 2 , the rate integrand Eq. (15) is weighted toward small q. In Fig. 8 we show the molecular form factors multiplied by F 2 DM ; the rescaled form factors are peaked more strongly towards low momenta, as expected. Because the inner peaks are kinematically forbidden for DM of all masses, but the tails of these peaks are also probed by all DM masses, this increases the magnitude of the peak-to-trough modulation amplitude to 70% for m χ = 2 MeV and remains as large as 30% for m χ = 100 MeV.

B. Target selection for daily modulation
The analysis of Sec. VI A suggests a strategy for designing target materials to obtain a large anisotropic response to electron scattering, and correspondingly large rate modulation, even in the limit of heavy DM. In this limit, the time-independent part of the argument of v − (Eq. (13)) is simply ∆E/q. Now consider a material with a form factor peaked at a momentum q * and with a lowest-lying excitation energy ∆E. To maximize the modulation, we look for a material for which q * and ∆E are related by q * ∆E vmax , where v max = v esc + v ⊕ is the maximum DM velocity attainable in the lab frame. Equivalently, the "effective velocity" characterizing the lowest-lying molecular transition, defined as v * = ∆E/q * , should be v * v max . In t-stilbene, the primary outer peaks have v * 200 km/s, which is a factor of a few too small to lead to the maximal rate, whereas the secondary inner peaks have v * 1200 km/s, and these peaks are always kinematically forbidden.
An ideal target for daily modulation would have either larger ∆E or a larger spatial extent (smaller q * ), so as to match v * v max for the primary peaks. To illustrate this, Fig. 9 shows the molecular form factor for the g → s 1 transition in t-stilbene but with the kinematically-allowed region defined by a transition energy ∆E = 8 eV, rather than the 4.2 eV in t-stilbene. Here, we have chosen a form factor F DM ∝ 1/q 2 , which weights the kinematically-forbidden inner peaks more, but the modulation is still driven by the forbidden region in q which has comparable radius to the outer peaks. As the forbidden region moves in q-space, the peak-totrough modulation amplitude can be as large as 20% for all m χ 20 MeV for the g → s 1 transition alone in this hypothetical material, almost a factor of 2 larger than the modulation amplitude for the equivalent transition in t-stilbene.
Taking a broader perspective, the anisotropic response of a condensed matter target to DM-electron scattering arises from an interplay of preferred scales q * set by the 14 FIG. 9. Same as Fig. 8 (bottom), but with a counterfactual transition energy of ∆E = 8 eV. This larger transition energy leads to a larger daily modulation for the g → s1 transition in the large mχ limit. To facilitate a comparison with Fig. 7, we also provide the modulation signal of the s1 transition with ∆E = 4.2 eV as the dotted gray line with the smaller amplitude. molecular size and a coincidence between the effective transition velocity v * and the maximum DM velocity in the lab frame. In the case of organic molecular solids, the conjugated π-electron system provides two length (momentum) scales given by the extent of the molecule along the molecular plane (q * 1.2 keV for t-stilbene), and the extent of a single 2p orbital (q * 6 keV), which sets both the carbon-carbon bond length and the extent of the out-of-plane π orbitals. The large hierarchy between these scales in large organic molecules means that the excitation dynamics in the plane are largely separated from excitations along the normal to this plane. Transitions along the normal direction will typically require larger imparted momenta than in the extended directions, and thus the form factor for the lowest transition will be peaked at q z = 0, with peaks in the x − y plane corresponding to the characteristic scales of the molecular (sub)structure. The kinematically-allowed regions which dominate the rate integral rotate in q-space over the course of the day, where the planar anisotropy (and to a lesser extent, the hexagonal structure of a benzene ring which breaks rotational symmetry to a discrete subgroup) gives the modulation for small m χ , and the anisotropy of the two displaced benzene rings contributes significantly to the residual modulation for large m χ . Having electronic transitions with v * slightly smaller than v max (as in our counterfactual example with ∆E = 8 eV in Fig. 9) will maximize the anisotropy for masses above threshold. That said, there is an inevitable tradeoff between the modulation amplitude and the total rate (consistent with the analysis of Ref. [60] for singlephonon production) because as v * approaches v max , the kinematically-allowed transitions rely more and more on the high-velocity tail of the DM velocity distribution.
From this perspective, we can understand why daily modulation amplitudes are typically small or nonexistent for electron scattering in conventional semiconduc-tor and noble liquid detectors. In noble liquids, the filled electron shells are spherically symmetric (ignoring small effects due to van der Waals attraction and dimerization between noble atoms), and thus the form factor will be isotropic and no daily modulation will occur. On the other hand, solid-state lattices have only discrete translational symmetries, which may be expected to lead to anisotropies like those due to the hexagonal structure of the benzene rings. However, the dominant lowenergy electronic transitions in conventional semiconductors with eV-scale gaps are due to delocalized valence electrons, which lead to a continuous energy spectrum and smooth form factors without a preferred momentum scale, at least for q smaller than the inverse lattice spacing ∼ 3 keV. 1 For larger q, scattering will probe core electron shells of single atoms at individual lattice sites, but these filled shells will be spherically symmetric and give isotropic form factors. That said, more exotic solid-state systems like Dirac materials, where a combination of a narrow gap (which permits small q) and an anisotropic linear dispersion ∆E ∼ v 2 x q 2 x + v 2 y q 2 y + v 2 z q 2 z with the v i bracketing v max , can have order-1 daily modulation [27,28] and a fairly large overall rate [16].
Importantly, q * is related to the characteristic size of the (sub)structure of the molecule as well as the symmetry of the transition, which determines whether the transition is dipole/quadrupole allowed. Meanwhile, the minimum ∆E is set by the HOMO/LUMO gap which is sensitive to the topology of the conjugated electron system, as well as the presence of functional groups which could donate or accept electronic density. For example, 1,2-diphenylacetylene is the acetylene-bridged analog of t-stilbene and presents a HOMO/LUMO gap of O(10%) larger than that of t-stilbene [89]. This implies that the two quantities are somewhat decoupled and can be independently tuned, at least for ∆E in the range 1 − 10 eV where the efficiency for scintillation photon detection is high. Furthermore, theoretical computations of the DM form factor as detailed in this paper can be verified by the complementary experimental probe known as electron energy-loss spectroscopy (EELS), which can be used to extract the generalized oscillator strength (i.e. dielectric function) of the molecular excitations in a given target; such a measurement automatically includes many-body effects [38] such as the ones parameterized by the PPP Hamiltonian as well as multi-electron excitations. Since t-stilbene and related chromophores have been identified as good scintillators for decades, their single-crystal synthesis is mature and can be scaled up to O(10 cm) crystals [90]. Therefore, it is entirely within the reach of existing methods and technology to implement O(10 kg) of organic crystal scintillation target since O(10 cm) single crystals of anthracene, t-stilbene and p-terphenyl are commercially available already, though ensuring the radiopurity of samples will be paramount to reduce backgrounds. In principle, this is no more challenging than obtaining radiopure liquid scintillator since crystals are readily grown from liquid stock.

C. Daily modulation from dark matter kinematic substructure
Another mechanism for exploring different scattering kinematics is supplied in the form of a cold, co-rotating stream of DM particles. The exemplary such subdistribution of DM particles is the putative Nyx stream [91]. This stream has velocity v Nyx (150, 0, 140) km/s [92] for components (v r , v θ , v φ ). In reality Nyx appears slightly anisotropic [91,92], but it has a relatively low spatially average velocity dispersion,σ Nyx 60 km/s. Given its inferred size, we can attribute to it a low escape velocity w esc = 150 km/s. To calculate the rate for the Nyx stream, v Nyx is subtracted from v ⊕ (40, 10, 230) km/s [92,93] in Eq. (13) and σ Nyx and w esc replace v 0 / √ 2 and v esc in Eq. (12), respectively.
Because of the smaller escape velocity w esc , activating the 4.2 eV transition requires larger momentum transfers, q O(10) keV. As a result the inner peaks of |f si | 2 at q 1.2 keV are kinematically inaccessible, and the peaks at q * 6 keV are only accessible at the high velocity tail of the distribution w ≈ w esc , even for m χ 100 MeV. Keeping the peaks of |f si | 2 at the edge of the kinematically-accessible region can induce a large modulation amplitude for a wide range of m χ , but at the price of significantly lowering the overall scattering rate. Because the Nyx fraction is 10% [93] of the local DM density, however, such modulation is unlikely to be a dramatic effect in any experiment, especially for the F DM ∝ 1/q 2 form factor.

VII. CONCLUSIONS
Among the many target materials proposed for DMelectron scattering, few have demonstrated the necessary anisotropic response to probe the daily modulation of DM, and none (to our knowledge) optimized for DM from the MeV to GeV scale. In this paper we have shown that organic crystals are a promising family of targets with excellent prospects for daily modulation, already at the same level as the expected annual modulation signal for the particular case of t-stilbene, and possibly larger if a compound with a suitable v * = ∆E/q * can be identified. In previous work [30] we have already demonstrated the efficacy of (liquid) organic scintillators in an experimental context, and we expect that many of the same design considerations will hold for solid-state scintillators.
The excellent overall sensitivity of t-stilbene -within a factor of a few of a comparable mass of silicon -combined with the additional handle of daily modulation, would make such a detector strongly complementary to the existing experimental program for Oscura [84] which uses silicon targets. In the event a positive signal is detected, daily modulation will be crucial for confirming a DM origin, and we have also derived a useful test statistic for determining the daily modulation significance in the presence of non-modulating backgrounds. We will explore design considerations for a concrete experimental implementation of a crystal organic scintillator detector in future work.
Beyond the particular case of t-stilbene, we have argued that aromatic organic crystals are a near-optimal compromise between overall rate, daily modulation, and scalability to large target masses, for DM of mass m χ 1 MeV. The building blocks of organic scintillators, the sp 2 -hybridized carbon orbital and its 2p z double bonding counterpart, are naturally anisotropic and support delocalized electronic states extended in the molecular plane, while the spatial extent of the 2p orbitals determines a preferred momentum scale. The weak intermolecular forces in organic crystals allow the electronic wavefunctions (and hence the form factors) to retain their molecular character rather than being entirely delocalized as in semiconductors. Furthermore, the discrete transitions at well-defined energies ∆E, combined with the sharplypeaked form factors, give v * which is close to optimal for t-stilbene, and may give larger modulation in compounds with slightly larger HOMO/LUMO gaps and therefore with a slightly larger v * . The combination of exciting features demonstrated by these results point to the ability to probe extremely well-motivated parameter space with plausible near-future technology, even in the presence of realistic but significant background rates. This indicates great potential for anisotropic organic scintillator detectors.
FIG. 10. Form factors |fs 3 | 2 and |fs 4 | 2 , shown in the qz = 0 plane as a function of (qx, qy). The s4 form factor has a roughly hexagonal structure, like the s1 transition, but stretched in the ±qy directions. The s3 form factor has approximately rectangular symmetry, stretched along the long axis (x =L) of the molecule. These two transitions provide the largest corrections to the scattering rate, but remain subdominant to the g → s1 transition even in the large mχ limit.
for the distribution of events in each bin. Defining for convenience, the mean µ 0 , variance σ 2 , skew and excess kurtosis of the Skellam distribution are [98] µ 0 = µ ∆ , so that in the large µ tot limit the distribution is approximately Gaussian. An exact version of the test statistic can be derived from the double-sided distribution, where the index (h) refers to the null or modulating hypotheses, (0) or (m). The difference between the test statistics, quantifies the significance of a signal and obeys a χ 2 distribution with two degrees of freedom. If µ 1,2 3, the Skellam distribution is well described by the Gaussian with σ 2 = µ tot , where the higher moments γ 1,2 become negligible. In this limit λ can also be approximated by an integral over a continuous variable, The significance of a measured ∆N can be easily expressed in terms of a number of standard deviations N σ by inverting the error function: As an example, we apply this result to a modulating signal, µ where f 2 is the integrated modulation fraction defined in Eq. (18), and the significance of a measurement of the number of events in each bin N 1 and N 2 can be assessed using Eq. (B9). In terms of R(t) from Eq. (15), µ s is where R is R(t) averaged over one sidereal day, and T exp is the total exposure time.
If the background rate were well understood, the total number of events (N tot = N 1 + N 2 ) could be compared to the prediction from the null hypothesis, N (0) tot = 2µ 0 b = µ tot as a way to discover or exclude particular DM models. Even without knowledge of the background, a small value for N tot can still be used to rule out those models which predict significantly more events than the measured N tot , but a larger N tot cannot be construed as a detection of DM without a better understanding of the background. However, the existence of a modulating signal provides an additional statistical handle on both discovery and exclusion.
Assuming that the background rate is unmodelled, a measurement of N tot supplies the best-fit values for µ in the null and modulating hypotheses, through All of the information about the signal, µ s , is extracted from the measured value of ∆N = N + − N − , which has expected value and we have explicitly specified that both the modulation fraction f and the signal strength parameter µ s depend on the DM mass m χ and cross section σ e . In assessing the capabilities of a directional detection experiment in the presence of daily modulation, we ask two questions: which DM models predict modulation signals that are large enough to be detected by the experiment? And, if the experiment measures a null result, which DM models are ruled out? The first question, which determines the discovery significance, can be posed in the context of ruling out the null hypothesis, where µ (0) ∆ = 0: in the Gaussian limit N tot 3. Models that satisfy N disc. σ ( ∆N ) > 3 or N disc. σ ( ∆N ) > 5 for this central value, for example, are likely to generate a modulation signal strong enough to claim a detection at the 3σ or 5σ level, respectively. The "Discovery" regions in Fig. 4 and Fig. 5 show N disc. σ (2f 2 µ s ) ≥ 3 using this central value.
To set exclusion limits, we ask which models are ruled out by a null result, ∆N ≈ 0. In this case, we use µ The exclusion curves of Fig. 4 and Fig. 5 show the models that would be excluded at the 90% confidence level from a null result ∆N = 0, using N excl. σ (0) > 1.65. In both examples above, we have taken ∆N to be equal to its expectation value under either the modulating hypothesis (∆N = 2µ s f 2 ) or the null hypothesis (∆N = 0), which has allowed us to quantify the significance of a measurement without reference to the log-likelihood ratio: both discovery and exclusion limits are given by This is simply because the "p value", λ, is equal to 1 at the central value of the double sided distribution Eq. (B5), and so the likelihood ratio is trivial. Because µ s is the expected counts per half of a day, we see that 2µ s is the exposure times the rate expected in Eq. (15). Thus, Eq. (B16) exactly recovers Eq. (19).
To assess a measurement away from the central value with more generality, it is better to use the log-likelihood L defined exactly in Eq. (B5) for each hypothesis. For small λ, the N σ (λ) defined in Eq. (B9) is approximated by This expression is particularly useful in the Gaussian limit, where N σ is given by Eq. (B9). To leading order in large N σ , where the test statistic ∆L compares a specific modulation model with f 2 (m χ , σ e )µ s (m χ , σ e ) to the null hypothesis, µ (0) b = 1 2 N tot . At the central values of the two distributions, ∆N = 2f 2 µ s and ∆N = 0, the test statistic takes values of ∆L = ∓4f 2 2 µ 2 s /N tot , respectively, and we recover exactly Eq. (B16).
For other values of ∆N , still in the Gaussian limit (µ s 3), the significance is found from the cumulative distribution function (CDF) of the χ 2 2 distribution, where γ is the lower incomplete Euler gamma function, and where k = 1 for the simple two-bin analysis. Generalizing to k statistically independent pairs of bins, the combined test statistic satisfies a χ 2 distribution with 2k degrees of freedom, and its CDF is given by In the limit of very few events, µ s 3, the CDF should be evaluated using the Skellam distribution instead, as it ceases to be approximately Gaussian for µ 1,2 < 1. However, as it is not possible to resolve an O(10%) modulation fraction with so few events, in this limit a more powerful constraint will come from using Poisson statistics on the total number of events.

Alternate Derivation with Poisson Statistics
The negative log-likelihood for a Poisson process is [ν k (θ) − n k + n k ln(n k /ν k (θ))] , (B23) where ν k is the expected number of events in a bin k, n k is the observed number in that bin, and θ are N θ parameters that determine ν k . The statistic L follows a χ 2 distribution of N bins − N θ degrees of freedom [99].
If we want to compare the hypothesis of a modulating signal versus the null hypothesis of a nonmodulating signal, we simply take the difference of their log-likelihoods. Since the total number of events in a day is fixed in the two scenarios, and only their distribution throughout the day is varying, we have where ν m k is the number of expected events assuming a modulating signal and ν 0 k is the number of events assuming a constant rate over the course of a day, and we have chosen the sign such that ∆L < 0 means that modulation is preferred. The distribution of values of |∆L| follows a χ 2 distribution with the number of degrees of freedom set by the difference between the number of parameters θ m and the number of parameters θ 0 . Eq. (B24) is an exact expression for the improvement in fit when allowing a modulating signal instead of a constant signal, appropriate for whichever event binning is most convenient. Because the total number of events is fixed, this is also the difference of the Kullback-Leibler divergences between these two hypotheses with the data.
Eq. (B23) and Eq. (B24) are appropriate for any choice of data binning, and can even be used for an unbinned analysis. For simplicity, and to provide an alternative derivation of the results in Sec. B 1, we calculate Eq. (B24) for the specific choice of two 12-hour bins per day, here labeled by + and −. The rates per bin are ν m ± =ν s (1 ± f 2 ) + ν b and ν 0 ± =ν s + ν b , where ν b is the expected background rate andν s is the average expected signal rate per bin.
Let us assume first that the true signal is not modulating: the number of observed counts in each bin in a given day is expected to be equal, such that n k+ = n k− = ν s + ν b . In this case, we have where in the second step we take the limit f 2ν s /(ν s + ν b ) 1. We now define N tot = d 2(ν s + ν b ) = 2N days (ν s + ν b ) to be the total number of events observed and, to make contact with the preceding section, we define µ s = dν s = N daysν s to be half of the total number of signal events expected to be observed over the entire experimental exposure. This gives The observed significance of a signal is therefore χ 2 = (2f 2 µ s ) 2 /N tot . Conversely, a limit at N σ significance on a modulating signal is possible when N excl. tot (2f 2 µ s /N σ ) 2 . As in the preceding section of this Appendix, because µ s is the expected counts per half of a day, the factor 2µ s is the exposure times the rate expected in Eq. (15). Thus, Eq. (B26) exactly recovers Eq. (19). Assuming on the other hand that the signal is modulating, the number of observed counts is no longer expected to be the same in the two bins in a given day. Instead, the counts will be related by n k± =ν s (1 ± f 2 ) + ν b . In this case, we have where the relative sign between Eq. (B25) and Eq. (B27) is reflective of the choice we made that ∆L < 0 means that a modulating signal is preferred. The magnitude of the significance is exactly the same as in the prior case, differing only in that the interpretation in this scenario is as discovery of a signal.