Dark Matter-Electron Scattering from Aromatic Organic Targets

Sub-GeV dark matter (DM) which interacts with electrons can excite electrons occupying molecular orbitals in a scattering event. In particular, aromatic compounds such as benzene or xylene have an electronic excitation energy of a few eV, making them sensitive to DM as light as a few MeV. These compounds are often used as solvents in organic scintillators, where the de-excitation process leads to a photon which propagates until it is absorbed and re-emitted by a dilute fluor. The fluor photoemission is not absorbed by the bulk, but is instead detected by a photon detector such as a photomultiplier tube. We develop the formalism for DM–electron scattering in aromatic organic molecules, calculate the expected rate in p-xylene, and apply this calculation to an existing measurement of the single photo-electron emission rate in a low-background EJ-301 scintillator cell. Despite the fact that this measurement was performed in a shallow underground laboratory under minimal overburden, the DM–electron scattering limits extracted from these data are already approaching leading constraints in the 3–100 MeV DM mass range. We discuss possible next steps in the evolution of this direct detection technique, in which scalable organic scintillators are used in solid or liquid crystal phases and in conjunction with semiconductor photodetectors to improve sensitivity through directional signal information and potentially lower dark rates.


I. INTRODUCTION
Dark matter (DM) can interact with electrons in a wide variety of systems, leading to a rich phenomenology of detector signatures and an active research program for development of new experiments . In particular, the formalism for DM-electron scattering in atoms [1] and solid-state systems [6,29] has been well-studied, but rather less attention has been devoted to DM-electron scattering in molecules. 1 In principle, molecules are promising detector candidates because covalent molecular excitation energies can be comparable to semiconductor band gaps, O(eV), allowing sensitivity to DM down to the MeV scale. Furthermore, achieving a large target mass (tens or hundreds of kg) of high-purity solvent such as benzene is somewhat easier than achieving a similar target mass of high-purity silicon, and when these molecules are used as solvents in a scintillating compound, the total background rate is quite competitive compared with the state-of-the-art in silicon achieved by SENSEI [16,20,23], CDMS-HVeV [21], and DAMIC at SNOLAB [24]. Finally, because the scintillation signal is decoupled from the primary DM-electron scattering event, any photodetector sensitive to the wavelength of scintillation light may be used to read out the signal. This is in marked contrast with many DM-electron scattering proposals to date, including silicon CCD's, where * carlosblanco2718@uchicago.edu † collar@uchicago.edu ‡ yfkahn@illinois.edu § blillard@illinois.edu 1 See [31], however, for detailed studies of DM absorption in molecules, and [32] for DM-nucleus scattering which excites rotational energy levels.
the target material itself serves as the detector.
In this paper we set the first limits on DM-electron scattering from an organic scintillator target, specifically EJ-301 [33], a ternary scintillator composed of 95% pxylene (1,4 dimethylbenzene) by mass, the rest being naphthalene and a proprietary fluor. We develop the theoretical formalism for DM-electron scattering in aromatic molecules (i.e. conjugated π-electron systems), exploiting the fact that semianalytic parametrizations of the electronic states can be found using a linear combination of atomic orbitals (LCAO) model. We then apply this formalism to the specific case of EJ-301, using experimental data from [34], where a low-background 1.3 kg scintillator cell was operated in a dedicated shield under minimal overburden. The method described in [34] allows the dark count rate of the photomultiplier tube (PMT) to be subtracted from the readout, resulting in a residual single photo-electron (SPE) rate of 3.8 Hz. As this background subtraction is only statistical in nature, and not event-by-event, it cannot be used to claim discovery, but nonetheless the residual rate is low enough to allow us to set an upper limit on the DM-electron cross section of about 10 −34 cm 2 for 10 MeV DM scattering through a heavy mediator, within an order of magnitude of the world-leading limits set by DAMIC in this mass range. We note that our work is complementary to previous work on DM-electron scattering in scintillating targets [10], which focused on solid-state systems as opposed to molecular solvents. This paper is organized as follows. In Section II, we define our molecular model for aromatic compounds and determine the relevant electronic wavefunctions and excitation energies. In Section III, we use these wavefunctions to compute the expected event rate for a given target mass. We present our results in Section IV. In FIG. 1. The chemical structure of benzene (left) and p-xylene (right). Following the convention common in organic chemistry, vertices are taken to be carbon atoms, single lines are carbon-carbon single bonds, and double lines are carboncarbon double bonds. The additional horizontal lines in pxyelene represent two CH3 methyl groups. In this case only one of the resonant Kekulé forms are shown for each molecule, but the rings should be understood to be entirely conjugated π-systems. Section V, we discuss potential next steps in the development of organic scintillator targets, and we conclude in Section VI with some thoughts on how to scale this setup to achieve superior limits. Some technical details of the wavefunctions and molecular form factors can be found in Appendix A.

II. MOLECULAR ORBITAL MODEL
Dark matter can produce a detectable signal in a liquid scintillator by exciting an electron from the ground state into one of several low-lying unoccupied molecular orbitals, which fluoresce upon de-excitation. The emission lines are broadened by roto-vibrational energy sublevels, thermal motion, and solvent effects. Therefore, the emission spectra of a liquid scintillator will be a continuum with peaks corresponding to the electronic transitions to be discussed in this section. The emission spectrum for EJ-301 is presented in ref. [35].
A prediction for the scattering rate requires knowledge of the momentum space wavefunctions of the bound state electrons in the molecule. In this section we find analytic expressions for these wavefunctions using a linear combination of atomic orbitals (LCAO). We then include configurational interactions (i.e. electron-electron repulsion) to obtain experimentally accurate descriptions of these electronic states.
In EJ-301, the primary target which initiates the scintillation process is p-xylene, a substituted benzene derivative with a chemical structure shown in Fig. 1 (right). In molecules such as benzene, the ring contains alternant double bonds. The electrons which occupy the π-orbitals in these bonds are essentially delocalized and are free to move along the so-called aromatic ring made up of the conjugated π-bonds. 2 In six-membered aro-2 Single covalent bonds between carbon are formed when two electrons occupy the σ-bonding orbital between two sp 3 -hybridized carbon atoms. The σ-orbital lies along the axis between the nuclei. In a carbon double bond, the carbons atoms are in an sp 2 -hybridized state where both the σ-bonding orbital and the matic rings there are six electrons associated with the conjugated π-electron system. The electronic transitions responsible for scintillation are those in which electrons in the conjugated π-bonds of the aromatic ring are excited into higher energy molecular orbitals. Thus, we will restrict our characterization to the π-conjugated system whose Hückel molecular orbitals are linear combinations of six 2p z atomic orbitals, one from each carbon in the ring. Each such carbon contributes one electron into the system. Following the wellknown LCAO method [36], the Hückel molecular orbitals (HMO's), ψ i , are given by where c i are the coefficients to be determined and φ 2pz (r−R i ) are atomic orbitals of the Slater type and R i are the equilibrium locations of the carbon nuclei. The 2p z Slater atomic orbital is given by where a 0 is the Bohr radius, and Z eff = 3.15 is the effective nuclear charge of the carbon 2p z orbital [37]. The HMOs can be determined by diagonalizing the six-by-six Hamiltonian, which reduces to solving the following linear system, where H lm = φ l |H core |φ m are the Hamiltonian matrix elements. We follow a type of Hückel model in which only nearest-neighbor interactions are considered. The Hamiltonian is then comprised of two types of matrix elements, one diagonal and one off-diagonal energy. The diagonal elements, i.e. the onsite energy, is an empirical quantity which varies by elemental atom. The off-diagonal resonance integral is a measure of the nearest-neighbor nuclear interactions. The values for the coefficients are given in Appendix A 1.
Following the notation of [36], we label the six HMOs in order of increasing energy as Ψ 2 , Ψ 1 , Ψ 1 , Ψ −1 , Ψ −1 , and Ψ −2 , where in the ground state of the molecule the first three are occupied and the last three are unoccupied (see Fig. 1). The minimum electronic excitation energy is from the highest unoccupied molecular orbital (HOMO) to the lowest unoccupied molecular orbital (LUMO), analogous to the valence-conduction gap in a semiconductor. The coefficients for these Hückel π-orbital are occupied. The π-orbital is lobed above and below the molecular plane, and is less tightly bound than the σ orbital. eigenstates are given in Appendix A 1. In this construction, (Ψ 1 , Ψ 1 ) and (Ψ −1 , Ψ −1 ) are degenerate pairs. However, since the the multi-electron wavefunction must be anti-symmetric, it is given by linear combinations of Slater determinants, i.e. the normalized and antisymmetrized products of six HMO's. The ground state for benzene is then given by: where |Ψ 1 , ..., Ψ n | is the antisymmetrized product of the HMOs Ψ 1 , ..., Ψ n , and Ψ is the opposite spin state as Ψ. In these products, the order specifies the identical electron indexing which implies |Ψ Following the method of Pariser, Pople, and Parr [37][38][39][40], we now include the electron repulsion term in the Hamiltonian, which becomes: where H core is the core Hamiltonian which gave us the HMOs and the second term is the two-electron repulsion operator responsible for the configurational energy of the electrons. In order to obtain the energies of each transition, the electron repulsion integrals must be calculated using the many-body wavefunctions [37,41]. Here, we adopt a semi-empirical model in which the first singlet excitation energies are tuned to the experimental values for p-xylene [42], and the second and third singlet excitation energies are calculated from the parameters derived in the literature [37]. Note that since the sensitivity reach of the scintillator is dominated by the lowest-lying excitations it is not particularly sensitive to large uncertainties in ∆E s2 and ∆E s3 . The molecular orbitals for the first singly excited singlet states are given by the following [36,43], where are single electron singlet excitations with respect to the ground state. The second and third singly excited singlet excitation wavefunctions and energies are given in Appendix A 2. Following the notation of [36,38,42], the ψ s1 1 and ψ s1 2 orbitals transform in the B 1u and B 2u representations of the point symmetry group, while the degenerate ψ s1 3 and ψ s1 4 transform as E 1u . The rate of dark matter-electron scattering depends on the molecular form factor, which specifies the probability of transferring momentum q to the molecule while exciting an electron from the initial state ψ i to the final state ψ f (hereψ(k) are the momentum space wavefunctions). Since the molecular form factor is an inner product of single-electron operators, the molecular form factors between the ground state and the excited states is evaluated via allowing form factors over many-body wavefunctions, ψ, to be computed in terms of the single electron HMO's, Ψ.
Since the methyl substituents in p-xylene are σbonded, they do not affect the π-electron system to first order in our model, and so the predicted wavefunctions for the conjugated π-electrons in p-xylene are equivalent to those derived for benzene. We account for the presence of p-xylene in the scintillator by tuning the excitation energies to those listed in Eq. (6). Note that the procedure described in this section for obtaining the many-body wavefunctions of electrons in benzene can be used to describe any molecule with a conjugated π-electron system, including the aromatic compounds used in other types of scintillators.

III. RATE CALCULATION
Dark matter with sufficient kinetic energy can induce a transition from the ground state to one of the excited bound states of Eq. (6) (or Appendix A 2) with a scattering rate determined in part by the momentum space wavefunctions of the initial and final electronic states. Conveniently, the details of molecular physics factorize into the form factor f ij (q) of Eq. (9), so that the scattering cross section (σv) i→j for any of the possible transitions reduces to [6] where M free is the amplitude for dark matter scattering on free electrons; q is the momentum transferred from the dark matter; m χ and m e are the dark matter and electron masses; i and j label the initial and final states, with a difference in energy of ∆E ij ; and is the energy transferred from the DM to the p-xylene molecule for an incoming DM velocity v (we have approximated the reduced mass of the DM-xylene system as just m χ which is appropriate for sub-GeV DM). It is conventional to separate the momentum dependence of M free into a model-dependent form factor F DM (q), written in terms of the reduced mass of the DM-electron system, µ χe , and the cross sectionσ e evaluated at a specific value of the momentum conventionally taken to be q = αm e where α is the fine-structure constant. F DM = 1 corresponds to a heavy particle mediating the DM-electron interaction, while F DM (q) = (αm e ) 2 /q 2 corresponds to a light mediator. A shift in the coordinate of a function r → r − R, as in Eq. (1), simply adds a phase e −ik · R to its Fourier transform. This property makes the LCAO model an especially powerful tool for obtaining approximate analytic expressions for the molecular form factor. In the case of benzene or xylene, the linear combination of atomic orbitals produces whereφ(k) is the Fourier transform of the 2p z Slater type atomic orbital Eq. (2), and where B ψ is a prefactor defined in terms of the coefficients c ψ i for each of the six HMO's Ψ i , the expressions for which are given in Appendix A 1 [44,45].
The total scattering rate R expected in the detector depends on details of the dark matter velocity distribution, g χ (v): where ρ χ is the local dark matter density, and N B is the number of p-xylene molecules in the scintillator. We approximate g χ (v) as a spherically symmetric speed distribution, and upon defining the expected rate of scintillation photons in the detector is where we have inserted an efficiency factor ξ, which is the product of the radiative quantum yield of the scintillator (that is, the probability of a primary excitation yielding a scintillation photon which escapes the liquid) and the quantum efficiency of the photodetector. A measurement of R puts an upper bound onσ e , after assuming a particular form for F DM (q) and integrating Eq. (19). Note that this rate does not include a contribution from electron ionization (as opposed to excitation to a bound state); while including ionization would only increase our rate estimates, the dynamics of free electrons in liquid scintillators and the corresponding scintillation rate estimate are much more involved than our simple treatment here, so to be conservative we neglect this contribution entirely.
Taking advantage of the fact that the six carbon atoms in the benzene and xylene rings are coplanar, every B ψ (k) can be written as a function only of some k x and k y , allowing the p z integral in Eq. (9) to be evaluated analytically. We perform the remaining integrals in Eq. (9) and Eq. (19) numerically using a Python implementation of the VEGAS adaptive Monte Carlo algorithm [46]. Upon expanding |f ij (q)| 2 as the product of two independent integrands, the rate R from Eq. (19) can be calculated relatively quickly from the resulting seven-dimensional integral. 3 In Fig. 3, we plot the angular integral of the absolute squares of the molecular form factors for each of the nine electronic transitions, i.e. dΩ q q 2 |f ij (q)| 2 . We note two important features: (a) the form factors peak at a momentum transfer, q ∼ 1/a 0 ≈ 2.5 keV, as expected for atomic systems, and (b) perhaps surprisingly, the lowest-energy 4.5 eV transition is strongly suppressed, such that the effective gap is closer to 5.6 eV. As we show in Appendix A 3, the form factors exhibit strong directionality, which is washed out under the assumption of a spherically-symmetric DM distribution but would be relevant for the true DM velocity distribution in the Earth frame, for which the DM "wind" exhibits diurnal and annual modulation. This suggests the possibility, which we discuss further in Sec. V, of using benezene as a directional DM detector: indeed, organic molecules including one or more benzene rings can exhibit a liquid crystal phase, and by applying suitable electric fields, nematic liquid crystals can have the constituent molecules aligned along a particular direction over large distances. We leave a dedicated analysis of the possibilities of such a liquid crystal scintillator for directional DM detection to future work.

IV. RESULTS
In a previously reported measurement [34], a 1.5 liter (1.3 kg) cell of EJ-301 scintillator, specially developed to minimize internal radioactive backgrounds, was operated in a dedicated shield under a modest overburden of 6.25 meters water-equivalent (m.w.e). By reducing the temperature of the PMT photocathode it was pos-sible to subtract the contribution from the PMT dark current to the SPE rate, generating in the process new limits on proton scattering by few-GeV DM candidates [34,47]. An irreducible SPE rate of 3.8 ± 0.1 Hz was obtained. Note that since this measurement was achieved over 14 days of data taking, comprising 4.6 × 10 6 events, the Poisson uncertainty of ∼ 2100 events = 0.0018 Hz is far below the 0.1 Hz uncertainty coming from other sources such as temperature control. Thus, we set limits by simply scaling the 3.8 Hz rate. This setup was recently upgraded to include an improved temperature control and discrimination against delayed PMT afterpulses [34]. The SPE background nevertheless remained unaltered, strongly suggesting that the present reach of the method is already limited by the environmental backgrounds associated to shallow underground operation.
To facilitate comparison with recent DM-electron scattering results, we take ρ χ = 0.3 GeV/cm 3 and evaluate the expected rate using Eq. (19) taking the speed distribution to be Maxwellian with a sharp cutoff at the escape velocity of the galaxy, where we take the typical velocity of the Earth to be v E = 232 km/s, local mean speed v 0 = 220 km/s, and escape velocity v esc = 544 km/s [48]. 4 The quantum efficiency of the PMT photocathode [50], integrated over the EJ-301 emission wavelengths [35] is 21.5%, and we calculate the radiative quantum yield of the scintillator to be 77% (see Appendix A 4), which gives ξ ≈ 0.166. As discussed in [34], the light collection efficiency of the EJ-301 cell is compatible with 100%. Fig. 4 shows the DM-electron scattering limit from EJ-301, conservatively assuming that 100% of the irreducible 3.8 Hz rate is due to DM interactions. Factoring out the 21.5% quantum efficiency of the PMT, the input-referred background rate is 17.7 Hz. Despite the large background, the limit for both F DM = 1 and F DM = α 2 m 2 e /q 2 is quite competitive with the recent DAMIC [51] and SENSEI [23] constraints, even exceeding the limits from XENON 10 [13] and XENON1T [25] below 5 MeV due to the lower threshold (5.6 eV for the lowest unsuppressed excited state in EJ-301 versus the 13 eV ionization energy of Xe) and surpassing the prototype CDMS-HVeV run [21] for m χ > 3 MeV. A modest improvement in the background rate, which could potentially be achieved either with additional overburden to reduce the cosmic rate or by using a photodetector with a lower dark rate, could set world-leading limits on . The solid black curve shows the limit derived from the 3.8 Hz residual background of 1.3 kg of EJ-301 scintillator operated in a shallow (6.25 m.w.e.) laboratory [34], and the dashed black curve shows the potential improvements from a background rate of 0.1 Hz. Shaded regions show existing exclusions from DAMIC [51], SENSEI [16,20,23], CDMS-HVeV [21], XENON10/100 [2], and XENON1T [25], rescaled as needed to a common DM density of ρχ = 0.3 GeV/cm 3 . In the bottom panel we show the stronger limit for XENON10 following the analysis of [13].
DM-electron scattering in the mass range 2 -7 MeV. We illustrate this in Fig. 4 by assuming a hypothetical background rate reduction to 0.1 Hz, giving the dashed black curve. For the silicon σe we use the 1e − example from [6]. With equal exposures, EJ-301 is comparable to silicon for FDM ∼ 1/q 2 for DM heavier than 10 MeV, and for FDM = 1 performs somewhat better for mχ 5 MeV.
dates since their electronic structures are well described by the formalism described herein [36], and their properties as scintillators are well studied as is the case for naphthalene and anthracene [52,53]. It should also be noted that as the excitation threshold decreases, one can expect the form factor to have the bulk of its support at lower momentum since the excitation energy and location of form factor peaks are both inversely propositional to the characteristic length between nodes in the wavefunctions, which increase for polyacenes with increasing numbers of rings. Both of these factors would increase sensitivity to lighter DM. Since the first transition in our idealized benzene model is suppressed by symmetry, one should expect this suppression to be lifted for larger and more complex (less symmetric) molecules, which would further lower the effective threshold and improve sensitivity to light DM.
The variety of aromatic molecules such as polyacenes, trans-stilbene, substituted benzenes, and oligophenyl chromophores, which are known to be good scintillators, provides a wealth of options in selecting possible detector targets. By running more than one experiment with different targets it may also be possible to discriminate backgrounds inherent in the main scintillator target since the expected rate for a given DM mass and cross section will be different for each molecule. Since these organic compounds are relatively inexpensive and the target is distinct from the photodetector, using several targets is only a marginal extension to using a single target.
Furthermore, the anisotropic nature of the crystal phases of these aromatic compounds can be used to look for directional modulation of the DM flux, further improving sensitivity to the DM-electron cross section and background rejection capabilities. In the solid state, large single crystals of polyacenes can be synthesized readily for pure compounds and binary/mixed scintillators [54][55][56]. These binary solid state targets maintain very high light yields while the crystal structure is known to produce anisotropic scintillation responses [57][58][59]. Aromatic liquid crystals based on anthracene core moieties, which can be aligned dynamically with electric fields while maintaining fluorescence quantum yield above 40% and threshold around 2.5 eV, can also be made in the laboratory [60]. Using this technology, one could imagine constructing a detector which tracks the DM wind by a continuous modulation of the electric field rather than a physical rotation of a crystal.

VI. CONCLUSIONS
In this paper, we have shown that organic scintillators are appealing targets for DM-electron scattering. The impressive reach, even for a non-optimized experiment, is due to two main factors. First, kilogram for kilogram, light DM scattering in aromatic hydrocarbons is as efficient as scattering in silicon for DM heavier than 10 MeV, even exceeding its sensitivity in the case of the F DM = 1 form factor. We illustrate this in Fig. 5, which compares the theoretical reach for 1 kg-yr of each material assuming a 100% quantum efficiency for signal detection (we still take a 77% radiative efficiency for EJ-301). We note that we cannot directly compare the form factors for pxylene and silicon because the former consists only of discrete transitions, while the latter has a continuous band structure, so instead we integrate over the form factor and the DM velocity distribution, which is effectively a calculation of the rate per unit mass. 5 Second, the scintillation process separates the primary scattering event (electron excitation) from the detected signal (a propagating scintillation photon), allowing for a separation of the target material from the detector. The scintillation photon, with energy of 2.9 eV, is in the near UV and can easily be transmitted across an interface (for example, an optical window) from the scintillator cell to any chosen photodetector. Other signals, such as electron/hole pairs and phonons, involve considerably more engineering at the interface in order to transmit the signal efficiently.
The only requirement is that the area of the photodetector be comparable to the area of the optical window in the scintillator cell, for maximum light collection efficiency.
These two factors immediately suggest a straightforward way to improve the reach of a liquid scintillator search. 6 By reading out the scintillation signal with a low-background silicon photodetector, either the phononbased sensor used in CDMS-HVeV or the skipper CCD used in SENSEI, we can take advantage of the best aspects of both setups: the large target mass (and low cost) of liquid scintillator, coupled to the low dark rate of the photosensor. Indeed, the reason our limits are competitive with the small-scale silicon experiments is that the effective background rate per kg of target material is very similar; the two-electron rate in SENSEI of 4.27 × 10 −5 /pixel/day with an active mass of 0.0947 g is equivalent to 4.7 Hz/kg, which is comparable to the 13.6 Hz/kg intrinsic background rate per unit mass of the EJ-301 (i.e. before accounting for the 21.5% efficiency of the PMT). The SENSEI sensitivity can obviously be improved by increasing the active mass, but assuming our background is not intrinsic to the scintillator itself, we could achieve similar limits by coupling the existing SENSEI detector to a kg scintillator cell, without the necessity of scaling up to a kg-scale mass of CCDs. 7 Further improvements would be possible by multiplexing several scintillator cells to achieve 10 kg of target mass, which, with zero background, would already be sensitive to a variety of thermal and non-thermal production mechanisms for sub-GeV DM [61]. In future work, we plan to measure the intrinsic background of low-background EJ-301 (processed with ion-exchange resins for uranium and thorium removal) with silicon photodetectors in the NEXUS facility in the MINOS cavern at Fermilab [62], which has 300 m.w.e. overburden and a 15 mK dilution refrigerator which allows the readout stage to be at cryogenic temperatures to reduce the dark rate.
In conclusion, we believe that aromatic organic scintillators are a promising addition to the panoply of novel condensed matter systems suitable for DM-electron scattering. By leveraging the advances in low-background photodetectors and reducing the intrinsic background in the scintillator as much as possible, we propose that 6 We are grateful to Noah Kurinsky for suggesting the improvements we propose here. 7 A potential concern with using silicon CCDs is that, since every scintillation photon will generate a single charge in a CCD pixel, high-energy background events producing many scintillation photons would appear as simultaneous single-electron events in multiple CCD pixels. This could be mistaken for a large singleelectron rate, if the CCD charge integration time is too long. However, this can be avoided with a sufficiently high CCD readout rate. In SENSEI, the continuous readout integration time is 20 ms for 800 samples with an RMS noise of 0.14 electrons [20]. In shielded, low-background conditions, a majority of CCD frames would not contain signals from high-energy backgrounds in the scintillator, avoiding this concern at the expense of a modest dead time. When the core hamiltonian, H core , is diagonalized, the LCAO coefficients in Eq. (3) are given by the following [36], (A1) where each c j is a collection of coefficients which gives an LCAO eigenstate. The normalization of these coefficients is such that each HMO is normalized to unity.
When these HMOs are transformed into momentum space, the phase prefactors in Eq. (14) are given by the following: 3 ak x 2 cos ak y 2 + cos(ak y ) where a = 0.14 nm is the bond length of the aromatic carbon-carbon bond in benzene and p-xylene.

Second and Third Singlet Excitations
The second singlet excitations, in which a single electron is excited across the second lowest energy gap, are given by the following: Notice that the mixing of configurational excited states results in an energy splitting around the quadruple degenerate second energy gap. It is found that these are the only linear combinations that give such a splitting. Finally, the third singlet excitation is given by

Form factor details
In this section we discuss the momentum dependence of the molecular form factors, including their angular profiles. Although the angular dependence of the scattering rate is washed out for a detector in the liquid phase, a crystalline scintillator would provide an opportunity to use directionality to discriminate a dark matter signal from the background.
The nine lowest-lying transitions described in Eqs. (6), (A3) and (A4) share a common feature: their form factors vanish when the imparted momentum q is orthogonal to the plane of the benzene ring. Adjusting the orientation of the molecule with respect to the dark matter wind can thus have a strong effect on the total scattering rate.
For concreteness, we use a coordinate system where the six carbon nuclei in the benzene ring are located in the z = 0 plane at: x + a sin j π 3 − π 6 ŷ (A5) for j = 1, 2, . . . 6, where a = 0.14 nm is the bond length of the aromatic ring in both benzene and p-xylene. All nine form factors exhibit the φ → φ + nπ 3 and φ → −φ symmetries of the benzene ring.
For m χ < 20 MeV the scattering rate is driven almost exclusively by the 5.6 eV and 6.4 eV transitions, and in Fig. 6 we show how their respective molecular form factors vary with the polar angle θ q for fixed azimuthal angle φ q = π 6 . The 5.6 eV form factor is maximized at θ q = π 2 and φ q = n π 6 for integer n, where the momentum transfer q is parallel to the displacement between neighboring carbon nuclei: While the 5.6 eV transition is driven mostly by inplane scattering with q ≈ 6 keV, the 6.4 eV transition is weighted towards lower momenta, peaked around q ≈ 3 keV for most values of θ q . The 6.4 eV transition is also responsive to a broad range of larger momenta 8 keV q 15 keV, but only if the polar angle is relatively steep, with θ q ≈ 13 • as shown in the second panel of Fig. 6. Although the 5.6 eV transition does exhibit some highermomentum response at this θ q ≈ 15 • angle, for a broad range of momenta centered around q ≈ 16 keV, it affects the scattering rate to a lesser degree.
Individual form factors also show strong directionality in the φ q direction, but often in complementary ways that tend to average out. For example, adding together the |f ij (q)| 2 form factors for the two degenerate singleelectron excitations that contribute to the 6.4 eV transition, we find a nearly rotationally invariant form factor. However, the single 5.6 eV transition retains its directionality, vanishing for φ q = n π 3 for integer n, and reaching sharp maxima at φ q = π 6 + n π 3 . Fig. 7 shows these two form factors along a particular conical slice at θ = π 2 , where the momentum q lies in the plane of the benzene ring.
As demonstrated in Fig. 3, the presence of the 4.5 eV transition is essentially irrelevant for DM-electron scattering, due to the strong suppression of its form factor. For completeness, it is worth commenting that its angular profile is complementary to that of the 5.6 eV transition; that is, for integer n its maxima occur at φ q = n π 3 , and its form factor vanishes along φ q = π 6 + n π 3 for all values of θ q .
If the dark matter mass m χ is large enough to probe the higher momentum behavior of the form factor, then the presence of the secondary peak around θ q ≈ 13 • for momenta q 8 keV has potentially beneficial implications for directional detection of dark matter. In addition to the primary response around θ ≈ 90 • , an observation of the relatively sharp secondary peak at θ ≈ 13 • could The azimuthal angle is fixed at φq = 30 • , which in our coordinate system is parallel to the edges of the hexagonal ring. Both form factors reach maxima at θ = π 2 where qz = 0, and vanish in the θ → 0 (q → qzẑ) limit. The 6.4 eV transition has an additional peak at higher momenta when 10 • θq 20 • . Neglecting the perturbation from the methyl groups in para-xylene, the form factors inherit the φ → φ + π FIG. 7. Azimuthal profiles of the 5.6 eV (top) and 6.4 eV (bottom) form factors |fij(q)| 2 , for momenta in the qz = 0 plane. While the 5.6 eV transition shows strong directionality in φ, vanishing at φ = n · 60 • for integer n, the sum of the two 6.4 eV transitions is approximately rotationally invariant.

Radiative efficiency calculation
The radiative quantum yield of the scintillator is given by the following expression [63], where S is the total scintillation efficiency, P is the primary excitation efficiency, C is the energy efficiency (taken to be 2/3 for organic scintillators), em is the energy of the emitted scintillation photon, ω is the energy gap of the excitation, Y is the experimentally measured light yield of the scintillator, and E abs is the total absorbed energy from an incoming particle (standardized to 1 MeV by light yield measurements). EJ-301 is measured by the manufacturer to have a light yield, Y 301 , which is 78% that of anthracene, Y anth = 1.74 × 10 4 photons/MeV [63], and an excitation energy, ω 301 of 4.5 eV. Note that using the lowest excitation in EJ-301 keeps our estimate conservative. Anthracene has an excitation energy of ω anth = 3.15 eV [52]. The excitation efficiency of organic molecules with aromatic rings is, P ≈ 2 3 F π , where F π is the fraction of π-electron in the molecule [63]. Here, we adopt values of the excitation efficiencies for EJ-301 and anthracene of P 301 =0.098 and P anth =0.01, respectively [63,64]. We can now calculate the radiative quantum yield of EJ-301, Q 301 as a function of the documented radiative quantum yield of anthracene, Q anth = 0.68 [52,63],

Molecular Form Factor: Analytic Calculation
The task of evaluating Eq. (19) with the molecular orbitals defined in Eq. (14) is simplified by the fact that the benzene ring lies in a plane, allowing the dp z integral in Eq. (9) to be completed by contour integration. In this section we list the result, as well as the outline for an analytic method capable of providing |f ij (q)| for an arbitrary arrangement of conjugated 2p z orbitals aligned with the sameẑ axis.
As described in Section III, the LCAO momentum space molecular orbitals can be factored into a commoñ φ(k), which is the Fourier transform of the 2p z atomic orbital, and an orbital-specific prefactor B ψ (k), which depends on the geometry of the molecule and the coefficients c ψ i . In terms of these two ingredients, the molecular form factor f ij (q) is described by where the newly defined I z contains the p z dependence: where we have absorbed the factors of a 0 and Z eff into q = (q x Z eff a −1 0 , q y Z eff a −1 0 , q z Z eff a −1 0 ) and likewise for p. Introducing as a convenient way to describe the locations of the poles in the complex plane, contour integration in dp z produces This analytic result significantly reduces the effort needed for the numeric integration, which can now be done over a smaller-dimensional volume. With additional work, the remaining dp x dp y integrals can be expressed instead in terms of generalized hypergeometric functions, based on the series expansion of I z in the azimuthal direction defined by the coordinate transformation p x = ρ cos φ, p y = ρ sin φ.
For greater generality, for the remainder of this section we take the locations of the nuclei R i and the coefficients c j ψ in the LCAO molecular orbitals to be generic: we require only that every R i lies in the z = 0 plane. As we did for p, we use cylindrical coordinates q z , q σ and q φ for the momentum transfer q, with q x = q σ cos q φ and q y = q σ sin q φ .
Given the coefficients Z k of the series expansion we proceed to derive a general expression for the molecular form factor in terms of the Kampé de Feriét hypergeometric function.
With this coordinate choice, f ij (q) has the generic form (A14) where for each pair (m, n) we define R mn = Z −1 eff a 0 |R m − R n | as the dimensionless distance between the mth and nth nuclei, taking the angular separation of the two nuclei with respect to the center of the ring to be φ mn . A benzene-like ring with uniform radius |R| satisfies Lastly we have defined ϕ ≡ φ − φ mn + π 2 (A16) so that Eq. (A14) reduces to the integral form of a Bessel function when I z is replaced by its series expansion. Defining the resulting expression for the dϕ integral can be simplified to after some manipulation of the indices in the series expansions. The appearance of the sinc(γπ) is simply to avoid double-counting the γ = 0 case.
The remaining dρ integration depends on the functional form of Z n . Rather than tackling the most generic case, whereφ(k) is an arbitrary momentum space atomic wavefunction expressed in terms of a Gegenbauer polyno-mial and spherical harmonics, we continue to specialize to the 2p z orbital and we evaluate Z n from Eq. (A12) using a Taylor series.
As f = ρ 2 + 1/4 is constant with respect to cos(φ − q φ ), the Z n can be determined from derivatives with respect to g = f 2 + 2ρq σ cos(φ − q φ ): Z n = (ρq σ ) n 1 g d dg n I z g→f .