Resonant absorption of bosonic dark matter in molecules

We propose a new class of bosonic dark matter (DM) detectors based on resonant absorption onto a gas of small polyatomic molecules. Bosonic DM acts on the molecules as a narrow-band perturbation, like an intense but weakly coupled laser. The excited molecules emit the absorbed energy into fluorescence photons that are picked up by sensitive photodetectors with low dark count rates. This setup is sensitive to any DM candidate that couples to electrons, photons, and nuclei, and may improve on current searches by several orders of magnitude in coupling for DM masses between 0.2 eV and 20 eV. This type of detector has excellent intrinsic energy resolution, along with several control variables---pressure, temperature, external electromagnetic fields, molecular species/isotopes---that allow for powerful background rejection methods as well as precision studies of a potential DM signal. The proposed experiment does not require usage of novel exotic materials or futuristic technologies, relying instead on the well-established field of molecular spectroscopy, and on recent advances in single-photon detection. Cooperative radiation effects, which arise due to the large spatial coherence of the nonrelativistic DM field in certain detector geometries, can tightly focus the DM-induced fluorescence photons in a direction that depends on the DM's velocity, possibly permitting a detailed reconstruction of the full 3D velocity distribution in our Galactic neighborhood, as well as further background rejection.

Dark matter (DM), a form of nonrelativistic matter that amounts to 25% of the energy budget of the universe but does not appear to emit light, is by now the conservative option to explain a wealth of astrophysical and cosmological data that can otherwise not be accommodated for with the known interactions and particles in the Standard Model (SM). The motion of stars in galaxies, the velocity dispersion of galaxies in clusters, gravitational lensing by galaxy clusters, temperature anisotropies in the cosmic microwave background, baryon acoustic oscillation measurements, and early-universe structure formation; all point to a new form of matter that is largely inert save for its gravitational interactions. Many questions remain unanswered: What are the properties-mass, spin, parity-of the dark matter particle(s)? What are its nongravitational interactions, if any? How is it produced?
While there are many possible answers to the first two questions, the number of dark matter candidates dwindles once you focus on the ones with a realistic production mechanism. One such great DM candidate is the socalled WIMP, a type of particle which may be produced with the correct relic abundance in the early Universe through the thermal freeze-out mechanism, provided it has a mass and interaction strength close to the electroweak scale. Searches for WIMPs are still in full swing, but previous iterations of both direct and indirect detection experiments have come up empty, ruling out the simplest implementations of this paradigm.
Ultralight, weakly interacting bosons constitute another large category of DM candidates with a natural production mechanism. Light spin-0 particles generically appear as relics of inflation through field misalignment, while massive spin-1 particles acquire an abundance set by quantum fluctuations during the last inflationary efold, among other possible production channels. As long as they are sufficiently weakly coupled such that they never reach thermal equilibrium throughout the cosmological evolution, bosons as light as 10 −21 eV can be cold DM, i.e. behave as an inert, pressureless, nonrelativistic fluid.
These bosons can arise in many theories beyond the Standard Model. The most famous example is the QCD axion, a light remnant in a class of theories that can explain the smallness of the neutron's electric dipole moment. Topological complexity in string compactifications naturally gives rise to a plenitude of bosonic states, such as axions and spin-1 fields, also sometimes known as "dark" or "hidden" photons. These states may easily have extremely weak couplings to the SM, as well as very small masses. Scalar fields associated with the shape and size of extra dimensions, as well as those that determine fundamental constants in our vacuum, often called moduli, can also couple very weakly and be extremely light. All of these states are associated with the same ingredients that give rise to the string landscape.
When bosons lighter than 15 eV make up a significant fraction of the local DM energy density, their number density is so large that there are many of them per de Broglie wavelength volume. When that happens, their superposition can be described as a classical field oscillating at a frequency set by the mass, and a coherence time determined by the inverse energy spread, roughly 10 6 periods of oscillation. This field also exhibits macroscopic spatial coherence on a length scale of order its deBroglie wavelength, 10 3 times larger than its Compton wavelength. The amplitude of the field oscillation is proportional to √ ρ DM , where ρ DM is the local DM density.
In this work, we take advantage of this behavior to propose a novel class of DM detectors. We describe how DM can act as a laser that resonantly excites transitions in molecules when its mass closely matches the transition energy, thus utilizing the DM's temporal coherence. Resonant absorption of DM can excite molecules to a higher-energy state that is otherwise not thermally occupied. This excited internal molecular state decays via emission of a photon, which eventually impinges onto a sensitive photodetector. Our techniques are applicable to DM masses between 0.2 eV and 20 eV, and can probe a variety of DM candidates, including axions, dark photons, and moduli. Two experimental configurations are shown in fig. 4. Molecular gas is placed in a container capable of supporting moderately high pressures. In the "bulk" configuration, a fraction of the container walls are instrumented with large-area photodetectors, and the rest of the container boundary is lined with an optically reflective coating to retain the isotropic fluorescence. The second, "stack" container exploits the spatial coherence of DM to focus the fluorescence onto a much smaller photodetector.
As we will show, the proposed setups have great intrinsic energy resolution and other advantages which allow for efficient background rejection and signal discrimination. In sec. II, we review the dynamics of a two-level system under influence of a nonrelativistic wave, and the types of molecular states and transitions that can be excited by bosonic DM. Section III contains a more detailed description of our experimental setup and strategy, as well as a discussion of backgrounds and signal discrimination techniques. We provide estimates for the sensitivity of our setup to scalar, pseudoscalar, and vector DM candidates in sec. IV. Finally, we compare our methods to other proposals from the literature in sec. V. Appendix A shows that the calculations performed in the semi-classical approximation throughout this work give-on average-the correct results.

II. THEORETICAL OVERVIEW
In this section, we give an overview of the relevant resonant absorption theory (sec. II A), cooperative radiation effects (sec. II B), the energy eigenstates of diatomic molecules (sec. II C), and the types of the transitions between them given certain operator structures of the darkmatter interactions (sec. II D). Some of this material is not new; we provide it merely to set up notation and give a self-contained review. We use natural units with = c = k B = 1 throughout.

A. Resonant excitation of a two-level system
Consider first a single molecule with two internal energy eigenstates |0 and |1 of an unperturbed Hamiltonian H 0 under which they have a relative energy splitting ω 0 with |0 the lower-energy state. We may parametrize the most general state of this system as a two-parameter space in a superposition angle θ and a relative phase ϕ which jointly define the surface of the unit (Bloch) sphere. The unperturbed time-evolution equations of this system are simplyθ = 0 andφ = ω 0 , determining the interaction-picture states |0 = |0 and |1 = e −iω0t |1 . We furthermore assume the excited state has a radiative decay rate of γ 0 due to spontaneous emission of photons, which will drive θ to π. We are interested in the dynamics of this system influenced by a weakly perturbing DM wave, in particular resonant absorption, which may be effected by a harmonic interaction Hamiltonian δH(t) with a small but nonzero matrix element 1|δH(t)|0 = Ωe −iα1 cos(ωt + α 2 ). (2) with α 1 an extracted phase such that the Rabi frequency Ω is real, and α 2 an arbitrary phase. In sec. II D, we develop the necessary tools to calculate the dark matter's Rabi frequency Ω given a nonrelativistic interaction Hamiltonian with nucleons, electrons, and photons, for the molecular states classified in sec. II C. Here, we first derive the absorption rate given a certain Ω, allowing an estimate of the minimum detectable δΩ later in sec. III. Our treatment is valid for any two molecular levels, but the reader may keep in mind the analogous systems in NMR or ESR, of a spin-1/2 particle with gyromagnetic ratio γ in a static magnetic field B 0ẑ , as well as an oscillating one δB cos(ωt + α 2 ). There exists an exact mapping of the spin projection states |0 ∼ = | ↓ and |1 ∼ = | ↑ (with the interaction picture being the rotating frame), the Larmor precession frequency ω 0 = γB 0 and the Rabi frequency Ω = γ|δB ×ẑ|, and the position on the Bloch sphere denoting the direction of the spin expectation value. Often the classical interpretation of precessing spins gives a useful intuition about the dynamics of the system.
Suppose the system starts in the ground state |0 at t = 0, at which point the DM wave is turned on. In the interaction picture, where δH(t) contains a term proportional to |1 0| = e iω0t |1 0 |, the state vector |ψ(t) evolves as (for short times t > 0): In the second line, we used perturbation theory firstorder in δH(t), and discarded rapidly oscillating terms of frequency ω 0 + ω, which give small corrections upon integration. For times t 1/ max{|ω 0 − ω|, γ 0 }, notice that the system rotates into a partly excited state with θ(t) = π − Ωt, and that the phase of the DM wave is imprinted onto the molecule as ϕ(t) = π/2 + α 1 + α 2 + ωt.
Heuristically, the state vector keeps precessing into a more excited state at angular velocityθ = −Ω until t ∼ 1/|ω 0 − ω| when the DM wave and the molecule dephase relative to each other, or until t ∼ γ −1 0 , the 1/e-lifetime of the excited state, whichever is shorter, so that the maximum excitation probability is | 1|Ψ (t) | 2 ∼ Ω 2 / max{|ω 0 −ω|, γ 0 }. At late times, we therefore expect an equilibrium to be reached between DM absorption and photon emission, each at a rate of γ 0 | 1|Ψ (t) | 2 . This intuitive result is also borne out by the fully-quantized treatment in ref. [1], where the absorption rate Γ (1) abs (ω) for a single two-level system at late times t γ −1 0 was found to asymptote to: Γ (1) abs (ω) = γ 0 Ω 2 /γ 2 0 1 + 4(ω 0 − ω) 2 /γ 2 0 + O Ω 4 , e −γ0t/2 , (4) for an excitation field in a coherent state. In App. A, we show that over integration times of interest, eq. 4 gives the correct result even when the field mode is not in a coherent state, as well as when the average occupation number in the mode becomes so low such that the semiclassical approxation breaks down.
One would naively think that extending this result to N = np 0 V molecules, of number density n in a volume V , and probability p 0 of occupying the state |0 , would be trivial, with a total absorption rate of N Γ (1) abs (ω). However, this is not correct, because it ignores cooperative effects due to the fact that the molecules interact with common excitation and radiation fields. Firstly, the nonrelativistic DM wave is expected to be phase coherent over lengths of order 1/mv 0 with v 0 ∼ 10 −3 a measure of the DM's local velocity dispersion, so that roughly the same phase ϕ is imprinted on the molecules within a sphere of this radius. Secondly, the radiation from nearby phase-matched molecules generally interferes, constructively so for separations close to an integer number of photon wavelengths 2π/ω, and destructively so for separations equal to a half-integer number of wavelengths.
In the next subsection, sec. II B, we explore these cooperative effects in more detail. In general, we find that the total radiative width of N molecules is not simply N γ 0 , but Nrγ 0 , withr an "average cooperation number" that can easily be much larger than unity. This leads to the surprising result that the radiative width per molecule is and depends on the molecular density, on the DM's phase coherence structure, and on the container size and geometry. For example, for a rectangular container of thickness R z > 1/m and where the other two dimensions are very large compared to the DM's de Broglie wavelength, we derive in sec. II B that, for noninteracting molecules, with n 0 the number density at standard atmospheric conditions with a gas pressure P = 1 bar and temperature T = 273 K. We defined p 0 as the thermal occupation probability of the state |0 . Indeed, we findr 1 for all but the largest and most dilute containers (before taking into account decoherence). For a vertical stack of slabs, the average cooperation number can be even larger by an amountS, which in an optimal arrangement is equal to the number of slabs in the stack within a coherence length.
An electric dipole transition has a spontaneous emission rate γ 0 = |µ 1,0 | 2 ω 3 0 /3π, where µ 1,0 is the electric dipole matrix element between |0 and |1 . In a stack of thin slabs, this then implies a radiative width at larger of: γ rad 8np 0S |µ 1,0 | 2 3mR z (7) ≈ 5.5 × 10 10 Hz p 0S mR z n n 0 For a single slab,S = 1 by definition; for a large stack of slabs, it can be as large asS ∼ 1/v 0 (see sec. II B for details). Finally, we are led to a total steady-state absorption/emission rate of: We note that on-resonance, when |ω − ω 0 | γ rad , the absorption rate is density independent, while the width of the resonance scales linearly with number density, because γ rad ∝r ∝ n. Off-resonance, when |ω −ω 0 | γ rad , the rate of DM absorption and subsequent photon emission scales as ∝ n 2 , where it can be regarded as a coherent conversion of DM into photons.
Besides radiative broadening, there are several other sources of resonance broadening, all of which we have ignored up until now. In an initial DM search-rather than a verification procedure of a signal-the only other broadening mechanism of importance is collisional broadening, as it is the only other one enhanced at high number densities. We will treat both the DM signal's frequency width and Doppler broadening in sec. III D, and nonradiative quenching in sec. III B. Inhomogenous broadening is always negligible. The effects of molecular collisions can be complex, but are often well-modeled by the "impact approximation", whereby molecules are unperturbed between collisions, but any state gets dephased to a new random phase between 0 and 2π when it collides elastically with other molecules [2]. (Inelastic collisions, where internal energy is exchanged between the molecules, typically occur less frequently.) Under these assumptions, effects from collisions lead to a Lorentzian lineshape g col 0 (ω 0 , ω 0 ) = 2 , with the mean collision rate per molecule in terms of an average collision cross-section σ col , the temperature T , and the total molecular mass M mol [3,4]. The size of the elastic collision cross-section as defined above is typically somewhat larger than the geometric size of the molecule. In the numerical estimate of eq. 9, we used σ col ≈ 10 2Å 2 as a typical benchmark value, and also assumed M mol = 40m p and a number density n ≡ P/T = n 0 at standard atmospheric conditions, with a pressure of P = 1 bar and temperature T = 273 K.
In the NMR/ESR analogy, γ −1 col is the equivalent of the transverse spin relaxation time T 2 , the timescale over which a pair of spins dephase due to magnetic dipole interactions.
We can fold in the effects of collisions into the total absorption rate formula, via the convolution Γ abs (ω, with γ ≡ γ rad + 2γ col . We should also note that the dephasing caused by the collisions also modifies the radiative width to γ rad γ 0 + η coh (r − 1)γ 0 with η coh given by eq. 30. Even so, it is still roughly true that γ ∝ n, so that the absorption rate is density-independent on resonance, and scales as n 2 off-resonance. A molecular sample in the gas phase will be at finite temperature, so typically a number of states will have appreciable occupation probability p 0 . The initial thermal state usually has nonnegligible support over several rotational states (sometimes also excited vibrational states) of the molecule, each occupied with the Boltzmann probability with j, k labeling all possible molecular states, each with energy E j . One can view the total molecular population as several independent subpopulations j with total numbers N j = p 0,j nV . Furthermore, each subpopulation will be sensitive to dark-matter excitations at a set of transition energies {ω 0,j } (not just one ω 0 as assumed above) far above the temperature T . This set of transition energies will be different for each subpopulation. The dense discretuum of states makes gas-phase molecules an effective multimode resonant system.

B. Cooperative radiation
The increased per-molecule radiative width of eq. 5 is a phenomenon known as cooperative radiation, closely related to the phenomenon of superradiance first described at length in ref. [5]. Cooperative effects are particularly dramatic for the emission following absorption of bosonic dark matter rather than photons, due to the DM's nonrelativistic nature.
This type of cooperative emission can be understood in classical wave mechanics, as it is simply due to the principle of superposition. Going back to the analogy between the two-level quantum molecular system and a precessing spin in a magnetic field, if the spins in an NMR/ESR sample are precessing in phase, they produce an effective magnetic dipole moment N times larger than that of a single spin. The radiated power then scales as N 2 , leading to a radiation reaction force per spin that is linearly proportional to N , in direct analogy to the Abraham-Lorentz force on an oscillating charge. This damping force leads to a characteristic damping rateand thus emission rate-proportional to the density of spins, just as in eqs. 5 and 6.
In the rest of this section, we perform a semiclassical computation of the cooperation numberr, a dimensionless number which parametrizes to what extent the radiation of different molecules interferes. For independent emitters,r is defined to be unity. However, constructive (or destructive) interference means thatr = 1 + O(n) in general. We will find below thatr depends on the phase coherence structure of the excitation source, as well as on the container geometry. In addition, for certain container geometries, the coherent emission can be highly directional.
Suppose we have a container of volume V filled with N molecules at discrete positions x ∈ V at large and uniform density n = N/V . The molecules are coherently excited by dark matter, which causes each of them to radiate a field at a spacetime point (t, x) of the form: In reality, the molecules will of course emit photons; as a proxy, we use massless scalar radiation to illustrate the essential physics. Also for simplicity, we will ignore the kinetic-energy contributions to the frequency of the wave, and thus take ω = m. In eq. 12, q is an effective "charge" set by the coupling of the molecule to DM as well as to the radiation field. The phase α(v, x ) of the DM wave is random: for a single velocity component v of the whole DM ensemble it can be written as The DM wave is, however, an incoherent superposition of many waves of different v and randomly distributed initial phase α v . We expect the DM to be approximately virialized in our Galaxy, with the kinetic energy having a Maxwell-Boltzmann probability distribution, which translates into a velocity probability density where v 0 ≈ 235 km/s is a measure of the DM velocity dispersion [6][7][8]. We have boosted to the laboratory frame which moves through the Galactic rest frame at velocity v lab . This relative velocity v lab = v + v is a vector sum of the Sun's velocity v ≈ 220 km/s in the DM's rest frame, plus the Earth's orbital velocity v ≈ 30 km/s around the Sun [9,10]. More accurately, this distribution should have a velocity cutoff at the Galactic escape velocity v esc ≈ 550 km/s [11], but we shall ignore this complication, as it will not significantly affect the results below.
The expected field at (t, x) from a single molecule at x is zero trivially because the phase average of cos[α + . . . ] returns zero. Note that averaging over the initial DM phases α v automatically takes care of time averaging. The same must thus also be true for the total field from all molecules in the volume V : The total emitted energy density m 2 A 2 tot does not vanish, as the square of the total field has the expectation value: In the second line, we have separated independent radiation terms and "cooperative" radiation terms, the latter of which capture interference effects. To get to the third line, we replaced sums with integrals as x → V d 3 x n(x ), used eq. 12, and averaged over time/phase. We also took the far-field limit, with |x − x | L −x · x (which gets simplified further to just L in the denominator factors) for all x ∈ V , and wherex is a unit vector along the line-of-sight direction. In the last line, we wrote the answer as the square of the total field from a single molecule q 2 2(4π) 2 L 2 , times the number of molecules nV , times a "directional cooperation number": The correlation function g(x, x − y ) itself is defined as For the Boltzmann distribution of eq. 14, we can evaluate this exactly: One can attribute the rapidly oscillating cosine factor mostly to the photon path length difference between |x − x | and |x − y |, with a longer-range modulation proportional to v lab from the coherent phase imprinted on the molecules by DM. The much more smoothly varying gaussian exponential factor quantifies the phase correlation between the states of two molecules with a separation of |x −y | and a coherence length of 2/mv 0 , equal to 0.5 mm for m = 1 eV. A similar correlation function was found to describe signals of bosonic dark matter in pairs of detectors located at different points in spacetime [12]. A directional cooperation number r(x) greater than unity means that the radiation constructively interferes in the directionx. Take for example the simplest and also most spectacular case, in which all N molecules are located within a DM's Compton wavelength of each other, |x − y | 1/m, ∀x , y ∈ V . In that case, the correlation function g B (x, x − y ) is unity for all directions and all separations, so r(x) = N , with perfect constructive interference.
Another instructive example is that of a rectangular prism with sides of lengths R x , R y , and R z aligned with the x, y, and z axes, respectively. For f (v) as in eq. 14, we can calculate r(θ, ϕ), the cooperation number for radiation at an angle θ relative to the z-axis, and at an angle ϕ in the xy-plane, as: Above, we have defined the dimensionless lengthsR i ≡ mR i , and the quantities β x ≡ sin θ cos ϕ − v lab,x , β y ≡ sin θ sin ϕ−v lab,y , and β z ≡ cos θ−v lab,z . We also defined FIG. 1. Angular dependence of the logarithmic intensity emitted by a molecular gas at large cooperation numberr 1 in a slab-shaped container with dimensionsRx =Ry = 10 3 and Rz = 1, whose orientation is depicted by the blue prism. The intensity is exponentially peaked in the directions close to the normal vector of each of the faces, and is largest by far along the z-direction, perpendicular to the two largest faces of the container. For illustrative purposes, we picked v0 = 10 −2 and v lab = 0, while we also exaggerated the thickness of the slab.
the integral function The limiting cases quoted in the second line are correct up to O(v 0 ) fractional error or better. Around |β| = 0 and for largeR, the integral function falls of very steeply as I ∝ exp(−β 2 /v 2 0 ). For a prism smaller than a 1/m on all sides, we recover the previous result of direction-independent coherent emission with r(θ, ϕ) = 1 + nR x R y R z , up to O(1/N ) fractional corrections. More surprising effects occur when some of the dimensions become large. For a "slab" with e.g. R x , R y 1/mv 0 , we find to a good approximation: for θ 1. The full angular dependence of the radiation from a single slab is depicted in fig. 1. We find that the coherent piece of the emission is highly focused, with 84% of the coherent radiation contained in a cone of angular radius v 0 ≈ 0.78 × 10 −4 rad, and 99.5% within twice that opening angle. For random angles θ, ϕ ∼ 1 and R z 1/m, one would typically find a much smaller FIG. 2. Density plot of the logarithmic intensity as measured by a horizontal screen placed far above a thin slab in the xyplane, zoomed in on the directions close to the normal vector defined by θ = 0, which pierces the screen at the location of the dot shown. The angular radius of the main emission region is the expected DM velocity dispersion measure v0 ≈ 0.78 × 10 −4 rad. The vector shows the angular offset of the cone's center to the normal; its magnitude and direction are given by the projection v lab⊥ of the average DM velocity v lab onto the xy-plane in the lab frame.
cooperation number r(θ, ϕ) ∼ 1 + n/(m 6 R x R y R z ). The coherent radiation cone is very nearly perpendicular to the slab, with the center of the cone (θ c , ϕ c ) slightly offset from θ = 0 by an amount determined by the lab velocity through the DM halo: The dependence of the emission cone's size and center on the DM's velocity dispersion v 0 and average velocity v lab is illustrated in fig. 1. In this way, a measurement of the angular distribution of the coherent radiation constitutes a measurement of the DM velocity distribution, including precise directional information! For a thin slab, the coherent radiation can dominate the isotropic, incoherent fluorescence. The appropriate measure is the average cooperation number Whenr 1, most of the outgoing radiation will be coherent. More precisely, the fraction of the radiation that is coherent is 1 −r −1 . For the thin slab, we find: For the small v 0 under consideration, the second fraction has a maximum value of 0.72 atR z ≈ 2.3, with subleading local maxima of {0.21, 0.13, 0.091, . . . , ∼ 2/R z } atR z ≈ {9.2, 15.6, 21.9, . . . , }. Numerically, we then find that the coherent emission can easily dominate for a gas (in absence of decoherence due to e.g. collisions) in a slab-like container at standard atmospheric conditions and assuming (locally) optimal thickness R z , as we have previously shown in eq. 6. The fact that the cooperation number in eq. 27 is inversely proportional to the thickness at large R z indicates that cooperative radiation is a surface effect rather than a bulk effect. We expect the scaling ofr inversely proportional to the (smallest) linear size R of the container to hold for any simply-connected, convex container volume for R 1/mv 0 . We furthermore expect this scaling and the overall value ofr to be quite insensitive to the particular shape of the velocity distribution f (v).
For disjoint container volumes, it is clear that constructive interference of the radiation from any pair of molecules in separate containers can be ignored if the containers are separated by more than a coherence length 1/mv 0 . When the containers are closer together, it is more appropriate to replace the sums x ∈V with integrals V d 3 x n(x ) and a piecewise uniform density n(x ), such that Disjoint container volumes may thus be arranged such that they interfere with each other in a controlled way. An important example is that of a stack of N s identical thin slabs (from the previous example) placed at a regular intervals along the z-axis, their centers a distance L apart. When the spatial periodicity along the z-axis is an integer number of wavelengths L = k2π/m, k = 1, 2, . . . , then the radiation from each slab constructively interferes with that of another, at least in the z-direction and for slabs within a coherence length. The full calculation of the directional cooperation number for such an arrangement is quite tedious, but can be done numerically with eq. 28. Roughly, the outcome is that the cooperation number of the stack is enhanced relative to r(θ, ϕ) for the single slab calculated in eq. 27 by an amount S(θ, ϕ): At best, the average cooperation number can be enhanced by about a factor ofS ∼ 1/v 0 for a judiciously chosen slab separation L = 2π/m. The above analysis has so far ignored decoherence, as it has implicitly assumed that each identical molecular system has retained its internal-state phase imprinted on it by the DM wave. Elastic collisions between the molecules will scramble the molecular phases, destroying any constructive interference in the radiation and driving r → 1 if it occurs at too large a rate. We estimate the coherent part of the radiative rate to be suppressed by the coherence efficiency in the presence of collisional broadening of width γ col as in eq. 9. This efficiency factor quantifies the fraction of the radiation emitted coherently (i.e. by the second term in eq. 17) with angular dependence proportional to r(x) − 1. The rest will be emitted isotropically. To focus the majority of the total radiative emission, one needs η coh > 1/2 or (r − 1)γ 0 > γ 0 + 2γ col .

C. Energy eigenstates of diatomic molecules
In this section, we provide a short but self-contained review of the types of molecular states, focusing in particular on the simplest system-the diatomic molecule-for pedagogical reasons. The classification below is a prerequisite to illustrate the different types of transitions that can be induced by dark matter; readers familiar with this material can skip ahead to sec. II D. Our abridged treatment below is almost entirely based on ref. [13].
The full Hamiltonian of a neutral, diatomic molecule in absence of external fields is: Here M 1 , M 2 are the two nuclear masses and the electron mass, α the fine structure constant, Z 1 and Z 2 the nuclear charges, −i∇ N (R N ) and −i∇ en (r en ) the momentum (position) operators for the N th nucleus and nth electron. We neglect spin-orbit coupling and relativistic corrections for the rest of this discussion; the results in this subsection and sec. II D are only valid insofar as Z N 1/α. Due to the large mass splitting of the nuclear masses and the electron mass-m e /M N ranges from 10 −3 to 10 −5 for light and heavy nuclei-the Hamiltonian H 0 is separable into "fast" electronic motion and "slow" nuclear motion (the Born-Oppenheimer approximation). Its internal energy (i.e. discounting translational motion) eigenstates |Ψ k with energy E k can be written as with the wavefunction |Ψ k factorized into one for electronic motion |χ el k , one for nuclear vibrational motion |ψ vib k , and one for nuclear rotational motion |Y rot k . The separability is manifested in the fact that transitions in the rotational and vibrational states leave the electronic state unaltered to a good approximation, and that the vibrational and rotational motion are largely factorized from each other. Different electronic states generally have different effective vibrational and rotational Hamiltonians, however.
Integrating out the electronic motion gives rise to a Hamiltonian of the form where we wrote the nuclear kinetic energies in terms of the center-of-mass momentum operator −i∇ Rcm conjugate to R cm ≡ (M 1 R 1 + M 2 R 2 )/(M 1 + M 2 ), and the relative momentum operator −i∇ R conjugate to the internuclear separation vector R ≡ R 2 − R 1 . We also defined the total molecular mass M mol ≡ M 1 + M 2 and the reduced nuclear mass M ≡ M 1 M 2 /(M 1 + M 2 ). The relative kinetic energy can be written out as: where J 2 is the molecule's angular momentum operator Parametrically, electronic energy splittings are of order ∆E el ∼ α 2 m e , and so is the binding energy U (∞) − U (R e ). The equilibrium radius R e is of order the Bohr radius a 0 = 1/(αm e ). Vibrational energy splittings turn out to be of order ∆E vib ∼ α 2 m e (m e /M ) 1/2 , while rotational energy splittings are yet lower at order ∆E rot ∼ α 2 m e (m e /M ), as we will see below. This hierarchy of scales is important, as it endows the molecules with a large set of absorption lines that is finely spaced and uniformly distributed in terms of transition energy, as we will explore further in sec. II D.

Rotational states
In anticipation of the hierarchy of vibrational and rotational energies, we can take R R e constant and study just the rotational spectrum ("rigid-rotor" approximation). At fixed R, the Hamiltonian simplifies to where the energy scale B e = 1/(2M R 2 e ) is half the inverse moment of inertia of the molecule. The energy eigenstates of this operator are just the spherical harmonics |JM ≡ |Y JM (θ, φ) with J = 0, 1, , 2, . . . and M = −J, − J + 1, . . . , + J which have the (2J + 1)-fold degenerate energies E rot J = B e J(J + 1). [These results hold true if the electronic level has zero spin S and zero orbital angular momentum Λ projected along the z-axis (see sec. II C 3 for definitions). For Λ = 0 but S = 0, each J level gets split into (up to) 2S + 1 sublevels, "Hund's case (b)". For Λ = 0, "Hund's case (a)", the J levels get split into two of opposite reflection symmetry σ v , and J has a lower bound of |Ω| ≡ |Λ + Σ|, where Σ is the component of the electron spin projected onto the molecular axis.] For example, the ground state J = 0 and the first excited state J = 1 are split by 2B e 1/(M R 2 e ) ∼ α 2 m e (m e /M ), after approximating R e ∼ a 0 = 1/(αm e ). The rotational constant B e typically varies by O(1) for different electronic states, as they stabilize the nuclei at different equilibrium separations. The variation of B e among different vibrational states is much weaker, and can be modeled by a correction factor α e : with v the vibrational quantum number (see below). The minus sign indicates that higher vibrational states typically have lower expectation values 1/R 2 . From naive dimensional analysis, α e /B e ∼ (m e /M ) 1/2 . Other corrections to the rigid-rotor approximation include centrifugal distortion, usually denoted by another correction term of the form D e J 2 (J + 1) 2 which is even smaller (D e /B e ∼ m e /M ) and will be ignored hereafter.

Vibrational states
Allowing for a dynamical internuclear separation R, we can integrate in the vibrational spectrum. Taking the energy eigenstate wavefunctions to be Rθφ|ψ vib Y rot = ψ v (R)Y JM (θ, φ)/R, we find that the radial wavefunctioñ ψ v (R) obeys the Schrödinger equation At leading order, we can ignore the R-dependence of the third term, the centrifugal potential, since it is typically a small correction to U (R). In that case, the vibrational energy eigenvalues and eigenfunctions become Jindependent. For a bound diatomic, the potential U (R) has a minimum at some R = R e , near which it can be approximated by a harmonic oscillator potential: which has eigenstates |v with energies for v = 0, 1, 2, . . . Anharmonicities in the vibrational potential can be modeled by the Morse potential: which has closed-form eigenfunctions with the exact eigenvalues: where ω 2 e ≡ 2a 2 D e /M and ω e x e ≡ ω 2 e /4D e . In other words, one recovers a pure harmonic oscillator for dissociation energy D e → ∞ while keeping the spring constant fixed at k e = 2a 2 D e .

Electronic states
The electronic structure and wavefunctions are more complicated than those of vibration-rotation modes. Nevertheless, the electronic levels can be classified according to their symmetry structure.
Because of the mass hierarchy between the electron and the nuclei, we can take the internuclear separation R = Rẑ as parametrically fixed to point along the z-axis.
(Note that we are thus temporarily adopting moleculefixed coordinates that are co-rotating with the molecule, rather than the space-fixed coordinates of the previous section.) One can find the electronic wavefunction and energy parametrically as a function of R, with the minimum of energy being T e at the equilibrium radius R e of the level under consideration. Knowing this parametric potential energy function U (R) also allows one to compute the vibration-rotation quantities like e.g. ω e and B e , which in general depend on the electronic level. In the Born-Oppenheimer approximation with factorizable motion of the electrons and the nuclei, we can consistently expand the energy eigenvalues in those of the electronic Hamiltonian with rovibrational fine structure for each electronic level (the total level labeled by k): Conventionally, the T e,k are measured relative to the minimum of the lowest electronic state, which is taken to have T e,k = 0. A diatomic molecule exhibits cylindrical symmetry, so the projection of the electronic angular momentum L z around the molecular axis is a good quantum number when the spin-orbit coupling is small. Hence we can classify electronic energy eigenstates according to their (molecule-fixed) φ dependence: L z |χ el = L z e ±iΛφ |χ el = ±Λe ±iΛφ |χ el = ±Λ|χ el States with Λ = 0, 1, 2, 3, . . . are conventionally denoted by the letters Σ, Π, ∆, Φ, . . . . Reflections around any plane through the molecular axis are also symmetries of the system, so eigenstates with Λ = 0 can be labeled by whether their wavefunctions are even (+) or odd (−) under such reflections, which are denoted by the operator σ v . In absence of spin-orbit coupling, total spin S and spin projection S z onto the molecular axis are also good symmetries, with quantum numbers S and Σ that can be half-integers. Finally, homonuclear diatomics (those with identical atoms) give rise to potentials for the electrons that are symmetric about the molecule-fixed inversion i, which takes (x, y, z) → (−x, −y, −z). Hence the electronic states can be categorized into those that are even (g) and odd (u) under i. For example, a homonuclear diatomic state with S = Λ = 0 that is even under σ v and i would be denoted by 1 Σ + g , whereas a heteronuclear diatomic state with spin S and orbital angular momentum Λ = 1 would be denoted by 2S+1 Π. Finally, we note that the electronic wavefunction must return to a sum of separated-atom wavefunctions in the limit R R e , so there exists a one-to-one correspondence between molecular orbitals and combinations of atomic orbitals. For the purposes of this discussion, we shall ignore complications due to spin-orbit coupling and nuclear spin symmetries; more details can be found in [13].
In figure 3, we depict a typical set of potential energy curves U (R) for five different electronic states in a hypothetical diatomic molecule. Insofar as the relative nuclear motion factorizes from the electronic motion, one can then compute the vibrational wavefunctions ψ v (R) in the potential well defined by U (R) for each electronic state separately. In the Born-Oppenheimer approximation, vibrational states can be excited within one electronic level. However, transitions between different electronic levels are generally accompanied by a change in nuclear motion-and thus a change in vibrational quantum number. The probability that, in an electronic transition from |χ el i → |χ el f , the vibrational state |v i in the initial electronic level changes to the vibrational state |v f of the final electronic level, is determined by the so-called Franck-Condon (FC) factor: This factor largely controls the relative absorption rates (and radiative emission probabilities) among the vibrational states. In the case of fig. 3 for example, if the electronic state is excited from χ el i = 1 Σ + g to the first-excited state χ el f = 1 Σ + u , the vibrational quantum number v i = 0 is much more likely to change to v f = 4 than to v f = 0 because of the larger vibrational wavefunction overlap (which can be read off visually from the wavefunctions in fig. 3). Tabulations of these FC factors have been computed and calibrated against measurement for a vast number of molecules; we will reference them where used in the text. In general, there are no selection rules con- 3. Electronic potential energy U (R) curves for a (hypothetical) homonuclear diatomic molecule as a function of internuclear separation R for five different electronic states, labeled by their quantum numbers, X-1 Σ + g , A-1 Σ + u , a-1 Σ + g , b-3 Σ + u , and B-1 Πu, respectively. Ticks in the potential well of each curve indicate vibrational energy levels, e.g. labeled from v = 0, ..., 6 in the fourth excited electronic state B-1 Πu. Exemplary vibrational wavefunctions ψv(R) are plotted in gray (with arbitrary vertical units) for a few vibrational states, including the lowest three of the ground electronic state. Rotational energy splittings are not shown. The ground state X-1 Σ + g can be excited to the first excited state a-1 Σ + g by a monopole operator, and to A-1 Σ + u and B-1 Πu by a dipole operator. A spin-dipole operator can also cause transitions to the spin-triplet state b-3 Σ + u , in addition to those caused by a regular dipole operator.
trolling changes in vibrational quantum number in electronic molecular transitions, greatly enhancing the number of potential absorption lines, and thus accounting for the much richer spectroscopy of diatomic and polyatomic gases than that of monoatomics. This increased complexity does not come at the cost of calculability nor spectral resolution for sufficiently small molecules. In table I, we list the defining spectral characteristics of several of the electronic energy levels used throughout the text.

D. Types of dark-matter induced transitions
In this section, we classify the types of transitions according to the operator structure of δH, and derive how they act on nuclear and electronic wavefunctions. For each type, we derive the corresponding selection rules for rotational, vibrational, and electronic transitions, and give estimates of the expected transition matrix elements in small diatomic molecules.  I. Spectroscopic properties of select diatomic molecules in the ground electronic state χ el = X, as well as in a few exemplary excited electronic states χ el = A, a, B, b, . . . taken from Ref. [14]. Energetic quantities such as the electronic excitation energy Te, vibrational energy splitting ωe and anharmonic correction ωexe, rotational constant Be, and vibrationrotation correction constant αe are expressed in the conventional inverse photon-equivalent-wavelength units of cm −1 ; conversion to units of energy is achieved via the substitution 1 cm −1 ↔ 1.23967 × 10 −4 eV. We also quote the equilibrium radius Re and boiling point temperature (for ground states) at standard atmospheric pressure P = 1 bar.

Monopole transitions
The simplest operator structure occurs when the perturbing Hamiltonian δH is spherically symmetric in all respects, in the sense that it acts trivially on the angularmomentum eigenstates of the entire molecule, the electronic configuration, and all spin degrees of freedom. In that case, δH must be diagonal in these angular eigenstate bases, so any such "monopole" transition obeys the selection rules: We expect these rules to be obeyed insofar as spin-orbit coupling can be neglected, except for the conservation of total angular momentum (∆J = 0), which is exact. It follows immediately that pure rotational transitions cannot be induced by a monopole operator. Vibrational transitions can occur, for example by the operators These act trivially on the angular wavefunction, but can excite an initial ground vibrational state |v i = 0 to a first-or second-excited vibrational state |v f via the matrix elements: In the limit of a pure harmonic oscillator (ω e x e → 0), these are the only nonzero, off-diagonal matrix elements connecting |v i = 0 to excited levels. In the presence of anharmonicities (ω e x e = 0), these selection rules are weakly broken, and e.g. (R − R e ) can also excite the ground vibrational state to v f ≥ 2, with the matrix ele-ments: up to O(ω e x e /ω e ) 1/2 -suppressed fractional corrections. In gross sensitivity estimates using the matrix elements of eqs. 47-52, we will often use the parametric estimates: typically correct up to a factor of 2, and "solving" for M and ω e x e in favor of the vibrational frequency ω e . Monopole vibrational transitions are not accompanied by a change in molecular rotation, and thus do not come with much rotational substructure in the possible transition energies (no dependence on B e ): The only dependence on the rotational quantum number enters via the vibration-rotation coupling α e , which in the Born-Oppenheimer approximation is suppressed by a factor of O(m e /M ) 3/2 relative to the vibrational splitting ω e . Monopole operators that can cause electronic transitions include Both the inversion i and the reflection σ v operations commute with these operators, so electronic transitions follow the selection rules in addition to the ones mentioned previously in eq. 45. Typical sizes for transition matrix elements in nonrelativistic molecules can be estimated via NDA as with R e the equilibrium radius of the initial electronic state |χ el i , usually an O(1) number times the Bohr radius a 0 = 1/αm e . The off-diagonal matrix elements of δH 0 IV and δH 0 V can be related as in eq. 58 because they are both component terms of H 0 in eq. 31, which itself acts diagonally on energy eigenstates by construction. Monopole transitions from the ground electronic state to an excited electronic state can occur at many different energies due to rovibrational substructure: We can see that even though ∆J = 0, there is a significant amount of rotational substructure because the rotational constants B e of two different electronic states generally differ by an O(1) fractional amount. As explained around eq. 44, there are no selection rules associated with changes in vibrational quantum number for electronic transitions, introducing a large amount of vibrational substructure, which also holds true for electronic transitions induced by other types of operators, to which we turn next.

Dipole transitions
Transitions caused by operators of the form are perhaps the most familiar, since they include the ordinary electric dipole transitions from photon absorption.
Here, the unit vectork denotes a unit vector in a spacefixed (as opposed to molecule-fixed) direction, like the direction of a vector field or the DM velocity, that acts trivially on the wavefunction. The effective directionk will in general vary in time and space, but only over scales of order the coherence time and length of the DM field, both much larger than the relevant temporal and spatial scales of the molecule for the energies under consideration.
The operator proportional to R = RR can induce pure rotational transitions through diagonal action of R = R e on the vibrational state but nontrivial action ofR on the rotational state, such that J f M f |k ·R|J i M i = 0. These well-known matrix elements have selection rules ∆J = ±1 as well as ∆M = ±1 fork ·ẑ = 0 and ∆M = 0 fork =ẑ. The rotational transitions which increase the internal energy are those with ∆J = +1, and have transition energy: These pure rotational transitions are too low in energy for the experimental setup under consideration in this work, because B e (for any molecule) is lower than the thermal energy T even at liquid nitrogen temperatures. However, they form an important part of the substructure for vibrational transitions. The operator R = RR can also act nontrivially on the vibrational state, e.
The transition amplitudes are products of the vibrational matrix element in eq. 47 and the pure rotational matrix elements. Over long time scales, the unit vectork will change direction, so we directionally average the rotational matrix elements ac-cording to: withR x ,R y ,R z unit vectors in the x, y, z directions, respectively. These average square matrix elements can be computed to be ( Higher-∆v are only weakly allowed (cfr. eqs. [50][51][52] and occur at energies that are roughly integer multiples of the first energy gap. To leading order in the Born-Oppenheimer approximation, the operator H 1 I does not excite electronic transitions.
The operator δH 1 II ∝ n r en primarily causes transitions between electronic states. The components of r en that point along the molecular axis are unaffected by L z rotations, are even under σ v reflections, and odd under the inversion i. Components of r en transverse to the molecular axis transform as Λ = 1 states under L z rotations, and are odd under σ v reflections and the inversion i. We thus find the selection rules: The last line follows from the trivial action on the spin coordinates, so both the spin multiplicity S and the projection of the spin onto the molecular axis Σ are unchanged in absence of spin-orbit coupling. For transitions that are allowed by these selection rules, we can estimate the matrix elements with NDA to be parametrically of order: Alternatively, for δH 1 II the above matrix elements can also be measured via absorption or emission intensities in electric dipole transitions.
As for any electronic transition, there are generally no selection rules for changes in vibrational quantum number. The total rotational quantum number J has the selection rule ∆J = ±1 for ∆Λ = 0 transitions, whereas both ∆J = ±1 and ∆J = 0 are allowed for ∆Λ = ±1 transitions (for most initial J i ). We refer the reader to [13] for more details about rotational fine structure in electronic transitions. The bottom line is that the extra rotational fine structure makes the discretuum of transition energies for electronic dipole transitions even more rich as that of electronic monopole transitions, with lines occurring at energies for transitions from the ground electronic state to one particular excited electronic state with excitation energy T e,f .

Spin-dipole transitions
Dark matter may also interact via less-familiar operators which exhibit a simultaneous coupling to both the spin and the momentum of nucleons or electrons: Despite the wholly different operator structure and besides their action onto the spin coordinates, they are in many ways similar to the dipole operators of eq. 59, with analogous selection rules on the orbital wavefunction in absence of spin-orbit coupling. For simplicity of discussion, we will take σ N to be the nuclear spin operator on only one of the nuclear spins (e.g. if the other one were spinless). For two spins, one would need to consider a larger spin Hilbert space, and use an interaction Hamiltonian of the form The first operator can cause pure rotational transitions, as well as vibrational transitions with rotational substructure. In spherical coordinates, we have Let us now see how this operator acts on a given rovibrational eigenstate. We have: where When the radial part of σ N · ∇ R acts diagonally on the vibrational state, we have that v|1/R|v = 1/R e and v|d/dR|v = 0. In this case, the angular and nuclearspin action can still cause pure rotational transitions along with nuclear spin transitions: Above we have assumed that the nucleus has a total spin S N equal to an integer or half-integer, and of which the component along the z-axis has the possible values Σ N = −S N , −S N + 1, . . . , +S N . The selection rules for these transitions are: We distinguish between the spin-flip (from σ N,x and σ N,y ) and spin-preserving (from σ N,z ) transitions, as the former can receive linear Zeeman energy shifts in an external magnetic field. They otherwise have the same rotational energy splittings and selection rules as in eq. 60. In a thermal state, the nuclear spins are in an unpolarized mixed state. The transition rate from an energy level with J = J i to one with J = J f is thus proportional to | J f |σ · ∇ R |J i | 2 avg , which can be found by averaging the square of the matrix element J f M f |σ · ∇ R |J i M i over the spin directionsσ, averaging over M i and sum-ming over the possible values of M f . One finds that The operator ∇ R can also induce vibrational transitions through its off-diagonal action. As can be seen from eq. 71, it acts on vibrational states |v through the operator terms d/dR and 1/R. The former has a matrix element connecting the ground state to the first-excited state: all other off-diagonal matrix elements from the ground state being suppressed by the anharmonicity of the vibrational potential. The second operator term, 1/R, has smaller vibrational matrix elements, as can be seen from and using eq. 47 to find that the dominant 1/R matrix element is smaller than that of eq. 74 by a factor of 1/M ω e R 2 e ∼ O(m e /M ) 1/2 . Therefore, the dominant vibrational transitions obey the same ∆v = ±1 selection rules as in eq. 62. The accompanying rotational matrix elements are the same as well, cfr. eq. 61, as long as one remembers that changes in the orbital angular momentum projected onto any particular axis are correlated with simultaneous changes in the nuclear spin projected onto the same axis as in eq. 73.
The operator σ e · ∇ e mostly excites electronic transitions, similar to δH 1 II of the previous section. The operator ∇ e has the same transformation properties under L z , i, and σ v as r e , so δH 1S II obeys the same selection rules as δH 1 II on changes in the orbital part of the electronic motion, which are already listed in eqs. 63-65. However, the nontrivial spin structure of δH 1S II means that it can induce both spin-preserving (∆S = ∆Σ = 0) transitionsthe only ones δH 1 II can excite-as well as spin-flip transitions (∆Σ = ±1 and ∆S = 0, ±1): Such spin-flip transitions are highly suppressed in small molecules for the usual dipole transitions, for which a ground state of e.g. 1 Σ can normally only be excited to higher spin-singlet states of symmetry 1 Σ and 1 Π, whereas these spin-dipole transitions can excite the same ground state to the spin singlets 1 Σ or the spin triplets 3 Π. For transitions that respect the selection rules of eqs. 76-79, we expect electronic matrix elements of size: with ω 0 = |E f −E i | the transition energy. The first equality follows by virtue of the identity ∇ en = −m e [H 0 , r en ] for the nonrelativistic H 0 from eq. 31. The matrix elements for the spin-dipole transitions can thus be inferred from those of the regular dipole transitions.

III. EXPERIMENTAL SETUP
We describe the general detector requirements and molecular container configurations in sec. III A, along with a detailed discussion of signal detection in sec. III B, background levels in sec. III C, and signal discrimination techniques in sec. III D. configurations. Molecular gas (depicted by light blue volumes) is pumped into containers capable of supporting pressures up to 10 bar. In the Bulk configuration, DM absorption events yield isotropic, single fluorescence photons, whose paths are indicated by thick red lines. Reflective coatings (shown as silver colored sheets) lining the container boundary retain the signal photons until they impinge onto a large-area photodetector, displayed as yellow tiles at the top. The Stack configuration features a pattern of alternating molecular density in the form of multiple slabs. This container geometry and the spatial coherence of DM can produce cooperatively emitted photons nearly perpendicular to the slabs, making it possible to focus them onto a tiny photodetector by a lens. The photon direction is sensitive to the DM velocity vector projected onto the plane defined by the slabs. For illustrative purposes, we do not show shielding, cooling, or electromagnetic field and pressure control systems in either setup. We also left out the reflective coatings on the front and top faces of the Bulk setup.

A. Configurations and search strategies
We envision two configuration types to detect bosonic DM in the mass range between 0.2 eV and 20 eV: one is a "bulk" detector volume, the other a layered set of slabs in a "stack" arrangement. We depict both configurations in fig. 4, and summarize their specifications in table II, which will be explained below. We consider a prototype "Phase I" and an optimistic, ultimate "Phase II" of both the Bulk and Stack configuration. We denote the four versions as BI, BII, SI, and SII.
The molecules are kept in the gaseous phase so that they can be regarded as approximately independent subsystems with a discretuum of energy levels, each with a resonant response to near-monochromatic excitations; intermolecular interactions in the liquid and solid phases tend to produce a non-resonant continuum in all but a few special cases. Dark matter waves of the right frequency-i.e. DM particles of the right mass-can excite molecules from their ground state(s) in thermal equilibrium to nonoccupied, higher-energy states, which in turn can fluoresce or cooperatively emit single photons. This radiation is to be read out by sensitive photodetectors, and serves as the signal in our setup.
Two primary considerations are key in a DM detector design based on resonant absorption onto molecules in the gas phase: radiative efficiency and frequency coverage. We will first explain the relevant physics for both of these requirements, and then how the Bulk and Stack detector designs address each of them.
Radiative efficiency-Ideally, every DM absorption quantum leads to a detectable fluorescence photon, as op-posed to heat or fluorescence photons that are difficult to detect (more on that in sec. III B). The dominant channel for conversion of internal energy of small polyatomic molecules to heat is via two-body collisions wherein the excitation quanta in electronic or vibrational state are converted to rotational and/or translational kinetic energy, a process called radiative quenching. (Other nonradiative pathways through which the absorbed energy can dissipate include vibronic decays in large polyatomic molecules, and molecular dissociation in weakly bound ones.) By analogy to the collisional broadening rate in eq. 9, we can thus parametrize the collisional quenching rate as: with a quenching cross-section σ quench independent of density. For many excited electronic states, the inelastic, quenching cross-section is typically an order of magnitude smaller than the elastic cross-section σ col [15]. Precise data on quenching cross-sections for excited electronic states is scarce, so we will take σ el quench ∼ 10Å 2 as a benchmark value. Vibrational quanta are less easily quenched in diatomic molecules, an effect theoretically understood in the context of Landau-Teller theory [15,16], with a quenching cross-section that depends strongly on relative molecular velocities and thus temperature. Ref. [17] found that a wide array of polyatomic systems obey the following empirical relation for the quenching rate γ vib quench of vibrational quanta, to a precision of 50%: with the definitions of the dimensionless quantitiesμ ≡ µ/m p ,ω e ≡ ω e /eV, andT = T /K. Above, µ is the reduced mass of the colliding pair of molecules, so µ = M mol /2 for a single molecular species of mass M mol . At sufficiently low temperature T , the quenching rate becomes exponentially slow, as indicated by the orange line in fig. 6, with about a factor 10 decrease in the quenching rate of carbon monoxide by lowering the temperature from 100 K to 77 K. (The long relaxation time of the molecular sample means environmental background rejection of cosmic rays and radioactivity may become harder in some cases, as we discuss in sec. III C.) We define the overall radiative efficiency η rad as the ratio of the fluorescence rate Γ rad to the DM absorption rate Γ abs from eq. 8 In general, η rad is a complicated function of the DM energy ω, as there are many possible target states that can be excited by dark matter, each having their own radiative properties. Given that Γ abs (ω) is highly peaked for ω near any one out of a set of transition energies {ω 0 }, we can take the overall radiative efficiency at ω to be that of the pair of states |0 and |1 with transition energy ω 0 closest to ω. The radiative efficiency around any such target state |1 can be written as: where γ 0 is the radiative width for the process |1 → |0 in vacuum, and γ i is the radiative width of |1 → |i with |i = |0 any intermediate state with energy below that of |1 . The fraction of absorption quanta that are coherently radiated is given by the second term in eq. 84, and is equal to: which is the generalization of eq. 30 in the presence of other radiative channels. For large cooperation numbers r−1 1 such that other radiative decays can be ignored, the condition for most of the absorbed to be coherently radiated is thus still (r − 1)γ 0 γ col . Likewise, the fraction of fluorescence radiation to the intermediate state |i is The sum over all of the intermediate-state branching ratios, i f i , can become close to unity when the intermediate radiative decays dominate over quenching at low number density, and when the radiative decay rate γ 0 back to the initial state is small because of selection rules, lowr, and/or combinatoric factors. Frequency coverage-Any one molecule can be regarded as a multimode resonator capable of absorbing DM particles with energies ω near any one of the set of transition energies {ω 0 }. It is important that most of the "frequency gaps" among the adjacent ω 0 can be efficiently covered, either via scanning the {ω 0 } by tuning some external variable like an external electromagnetic field, or by broadening each individual line by increasing the molecular number density.
The splittings among the possible transition energies are typically of order the rotational energy constants of the ground vibrational state. A diatomic molecule has rotational energies E rot = B e J(J + 1), J = 0, 1, 2, . . . in its ground state, where the magnitude of the rotational constant B e is determined by the inverse moment of inertia, i.e. the reduced mass M of the diatomic times its mean square separation R 2 R 2 e : For example, for dipole vibrational transitions, we expect a rotational fine structure with a splitting of 2B e , as derived in eq. 62 and depicted in the lower panel of fig. 5. Most other types of transitions have similar fine structure splittings at the same order or lower. Polyatomic molecules composed out of more than two atoms typically have multiple rotational constants (one around each axis, unless it is a linear molecule) that are smaller, and thus exhibit an even richer fine structure and lower degeneracy. For example, SF 6 has rotational constants of 1.1 × 10 −5 eV around all three axes.
The splittings in rotational transition energies can be scanned by tuning the magnitude of an external electric field via the second-order Stark effect, if the molecule exhibits a permanent electric dipole moment µ e (in molecule-fixed coordinates, not laboratory-fixed coordinates of course). In an electric field E, the energy eigenvalues E rot (J, M ) of the combined rotational and Stark Hamiltonian are becomes O(1). In the above numerical estimate, we took the electric field to be the breakdown electric field value (for large separations) of air at standard atmospheric conditions, and a dipole moment of 1 D ≈ 0.39ea 0 . Transitions involving an electronic spin flip can furthermore receive linear Zeeman corrections to their transition energy in an external magnetic field B, by an amount with g e the effective g-factor and µ B the Bohr magneton. For nuclear spin flips, the Zeeman shift is of less relevance due to the smallness of the nuclear gyromagnetic ratio. We thus conclude that it is possible to scan the gaps in the rotational fine structure of the absorption spectrum for molecules with small rotational splittings, which include heavy (large M ) and weakly bound (large R 2 ) diatomics, or moderately large polyatomics (large moments of inertia, usually around several axes).
Another important way to achieve contiguous frequency coverage is to broaden all of the transition lines {ω 0 } at high number densities n. The full width in ω at half-maximum absorption rate in eq. 8 is where we took T = 273 K and a molecular mass of M mol = 40m P in the collisional width 2γ col defined in eq. 9. Collisional broadening alone can already bridge the gap between typical rotational splittings, cfr. eq. 87, at moderate pressures of 10 bar. The radiative width γ rad = γ 0 +η coh (r −1)γ 0 + i γ i can increase the FWHM even more, especially ifr n/m 3 and γ 0 is an E1-allowed decay rate, cfr. eq. 93. Bulk configuration-The Bulk detector configuration comprises of a large convex volume, V = (0.3 m) 3 in Phase I, V = (2 m) 3 in Phase II, filled with a single molecular species at any time. A large-area photodetector array (PMT in Phase I, MKID array in Phase II) detects the signal photons. The rest of the container area is coated with a highly reflective layer, necessary to retain the fluorescence photons since they are emitted isotropically. Cooperative effects may be ignored for a bulk container, i.e. η coh (r − 1) 1, as we showed in sec. II B. Photons from the |1 → |0 fluorescence channel are at a danger of being reabsorbed elsewhere in the detector volume, and are thus more easily quenched after multiple re-absorptions and re-emissions (we discuss optical thickness issues in more detail in sec. III B). Optical thickness is a problem whenever the matrix element of the |0 ↔ |1 is E1-allowed and γ 0 is thus relatively large; on the other hand, if |0 ↔ |1 is E1-forbidden, then γ 0 can usually only dominate over the quenching rate from eq. 81 at very low number densities, potentially suppressing the absorption rate.
A Bulk detector configuration will therefore largely rely on radiative decays |1 → |i where |i is any intermediate state with negligible occupation probability in thermal equilibrium. These decays will produce fluorescence photons for which the medium is optically trans-parent. Hence, the detectable photon fluorescence rate Γ B rad for a Bulk configuration is where it is understood that f i from eq. 86 is to be evaluated for the target level |1 with ω 0 closest to ω, and that the sum is to be performed only for intermediate levels |i that produce fluorescence photons for which the molecular medium is optically thin. Frequency coverage for a DM search over a broad range of masses in a bulk detector can be achieved by performing several narrow searches, each with a different molecule. Covering gaps in the rotational fine structure of the absorption spectrum can be done in the two ways mentioned above: scanning at low density or broadening the lines at high number density. Within the scanning strategy, the optimal regime is to work at a low number density where γ quench thin i γ i , such that the fluorescence fraction is large, The resulting narrow linewidth γ means that many (on the order of γ/B e ) shots, each at a different set of {ω 0 }, are needed for contiguous frequency coverage. Each shot can thus last a small fraction of the total integration time t int for the molecule: t shot ∼ (γ/B e )t int . At high number densities where γ = 2γ col B e , only one shot is needed, but at the cost of a lower fluorescence ratio thin i f i thin i γ i /γ quench ∝ 1/n. At finite background rate levels, the scanning method typically provides better optimum sensitivity. In the limit of zero background, they yield roughly equivalent sensitivity, though the broadening method has the advantage that it can be used for all transitions in all molecules (including those without permanent electromagnetic moments), and is experimentally more straightforward.
Stack configuration-To take full advantage of the cooperative radiation effects analyzed in sec. II B, the detector volume should be composed out of a planar stack of slab-like containers made out of highly transparent material such as glass or silicon, as shown in the bottom diagram of fig. 4. The signal photons are emitted nearly perpendicularly to the planes, and could thus be focused onto a much smaller photosensitive area by a lens or reflecting mirror (or an array thereof). Any deviation of the emission direction from the normal vector is of order the DM velocity divided by the speed of light, so a highresolution photodetector such as an MKID could "image" the DM velocity distribution as in fig. 2 as a function of time.
Alternatively, an "artificial stack", without physical barriers between the slabs, may be effectively created by a standing electromagnetic wave pattern in a bulk container. At the antinodes of the standing waves, the ground state population could be depleted via resonant pumping to an intermediate level, while the molecules would be much less affected near the nodes of the standing waves. Even if no such suitable intermedi-ate level is available, the standing wave pattern could create a spatially-dependent quadratic Stark shift, moving the molecules off-resonance at the antinodes and onresonance at the nodes (or vice versa) at sufficiently small line widths. This artificial stack approach likely introduces additional complications; an evaluation of its feasibility is beyond the scope of this work.
The aim of the Stack design is to operate in a regime where the equality is roughly satisfied. This accomplishes three goals simultaneously: radiation focusing, high radiative efficiency (η coh ≈ 1), and potential contiguous frequency coverage. The coherent, focused fluorescence rate is: The FWHM of the DM absorption lines will be as in eq. 91 with the radiative width dominating if eq. 93 is satisfied. Molecules with an E1-allowed transition with transition dipole moment µ 1,0 = 1|µ e |0 in a stackedslab configuration have a cooperative radiative width as in eq. 7, so we have that the ratio controlling the validity of eq. 93 is for σ col = 10 2Å 2 , T = 273 K, and M mol = 40m p . The radiation focusing effect analyzed in sec. II B has several important advantages. It allows for a smaller photosensitive area, permitting the use of highly sensitive, cryogenic photodetecors such as a TES or MKID without a prohibitively high cost. In addition, it isolates the signal from the mostly isotropic environmental backgrounds (see sec. III C). The direction of the coherent photons depends on the DM velocity and dispersion as shown in sec. II B, yielding a spectacular intrinsic sensitivity to the direction of dark matter. Low-dark-count, 10-kilopixel MKIDs as large as (10 mm) 2 have already been built, and would have an intrinsic DM velocity resolution of 10 −5 (roughly given by the pixel size divided by the transverse size of slab) in the two dimensions parallel to the stack, i.e. at the sub-percent level fractionally given the expected DM velocity at 10 −3 the speed of light. Full 3D velocity information could be gleaned over long integration times because of Earth's rotation and orbit, and/or with two experimental setups.
Sensitivity estimates-We quantify the sensitivity of each detector version in terms of the smallest Rabi frequency δΩ it can detect at unity signal-to-noise ratio (SNR) after a shot time t shot . A photodetector with an intrinsic dark count rate DCR within its bandwidth ∆ω around an energy ω has a SNR = 1 sensitivity to detected photon counts at a rate of where Γ bckg is the true photon background rate impinging on the photodetector with detection efficiency η det , and in the same bandwidth ∆ω. If we define η γ as the probability that a radiated photon is detected by the photodetector, then the signal detection rate is The minimum detectable Rabi frequency δΩ is the Ω for which Γ det = δΓ det . When the absorption rate at a given energy ω is dominated by a single resonance with transition energy ω 0 closest to ω, we can write: In table II, we have compiled the typical on-resonance sensitivity in terms of δΩ for p 0 = 0.1, γ = 2γ col , σ col = 10 2Å 2 , and M mol = 40m P for each configuration and phase, assuming Γ bckg DCR + t −1 shot . Comparison with the signal Rabi frequencies Ω in table III allows one to estimate the signal-to-noise ratio for many DM candidates at a few benchmark frequencies. We shall see in sec. IV that the Phase I prototypes already explore new parameter space of some DM candidates, and that Phase II experiments will be capable of probing previously unexplored parameter space in all of the DM models and couplings listed in table III.
In figures 6 and 7, we illustrate many of the considerations in this section so far with a simple case study, namely the first-excited vibrational level in carbon monoxide. Figure 6 shows the inverse time scales of radiative and collisional dynamics as a function of pressure and temperature. We also show on the same axis the rotational energy constant B e , and the rate Γ RD to which radioactive decays can be held in the large BII setup with low-contaminant materials (before any active veto). Figure 7 shows in green the absorption rate Γ abs of a kinetically-mixed hidden photon DM particle with = 10 −12 (see sec. IV A for more details) for the 12 C 16 O isotope in the SI setup. The resulting coherent radiation rate of eq. 94 is plotted in red, about 2-3 orders of magnitude below. This is because CO has a rather small transition electric dipole moment, and thus η coh 1 even with the assumed stack parameters ofS/mR z . Molecules with stronger absorption strengths-such as e.g. HCl, CO 2 , and CH 4 -can achieve η coh much closer to 1.

B. Photodetection
In sec. II A, we calculated the absorption rate of the DM particles of energy ω given a Rabi frequency Ω, or equivalently the rate at which molecular states with transition energy ω 0 close to ω are excited. The signal detection consists of reading out the fluorescence photons spontaneously or coherently emitted from these excited levels, at a rate Γ rad = η rad Γ abs . However, the photon Dynamic rates of the two-level subsystem |v = 0, J = 1 ↔ |v = 1, J = 0 in carbon monoxide (CO), as a function of pressure P . Plotted are the collisional broadening rate γ col (green), the coherent emission rate in absence/presence of decoherence-(r − 1)γ0 in dashed red / η coh (r − 1)γ0 in solid red-and the nonradiative quenching rate γ vib quench (orange), for two temperatures T = 100 K (thick) and T = 77 K (thin). Also shown are the incoherent radiative rate γ0 (blue), and a benchmark natural radioactive decay rate ΓRD (dashed purple) before an active veto. The collisional width γ col comes within a factor of 3 (30) from the rotational energy Be (dashed black) at the boiling point pressure of 5.5 bar (0.6 bar) at T = 100 K (T = 77 K). This shows that collisional broadening is an effective mechanism for frequency scanning. Notice that because the collisional rate dominates every other dynamical timescale, it reduces the coherent emission rate so that only 1 every O(1000) DM particles absorbed will produce a photon. detection rate is typically lower than this radiation rate Γ rad , cfr. 97, by the overall detection efficiency factor We describe our estimates for η γ below, and show that it can be O(1) in our setup. We will finish this subsection with a brief summary of potential photodetectors and their specifications, including dark count rates. Four main loss mechanisms are responsible for η γ ≤ 1, each with their own efficiency factor η i as schematically indicated in eq. 99. They are due to, respectively, absorption onto the reflective walls of the container with probability 1 − η refl , finite transmission probability η trans of glass elements in the optical path, re-absorption and subsequent fluorescence quenching (if the gas is optically thick) with probability 1 − η thick , and the intrinsic quantum detection efficiency η det of the single-photon counter.
Rate of absorption Γ abs (solid green) and rate of focused, coherent fluorescence Γ S rad (solid red) as a function of DM energy ω m, in the Phase I version of the Stack configuration, filled with carbon monoxide (the 12 C 16 O isotopologue) at a pressure of 5 bar, temperature of 100 K, and "stack enhancement parameter"S/mRz ∼ 10, cfr. eq. 7. The case study shown assumes a kinetically-mixed photon with = 10 −12 , zoomed in on the energy range around CO's first vibrational absorption line |v = 0 → |v = 1 at ωe = 0.27 eV, with rotational fine structure of splitting 2Be clearly visible. This plot shows that the thermally occupied rotational levels of the ground state allows us to expand the frequency coverage of vibrational and electronic transitions. In dashed green, we also show the absorption rate for 12 C 18 O, displaying the isotope shift on the rovibrational structure.
In what follows, we discuss how each of the η i can be of order unity.
Reflection efficiency-Fluorescence photons may be absorbed by the (not perfectly) reflective walls of the gas container in the Bulk configuration. We envision that the container walls are coated with a material with high reflectance R, which depends on the wavelength, polarization, and incidence angle of the incoming photon. Appropriately averaging over polarizations and incident angles for an effective wavelengthdependent reflectanceR(ω), a photon can be reflected an expected 1/[1 −R(ω)] times before being absorbed onto the coating. The loss fraction of signal photons will be small η refl ≈ 1, as long as A det /[V 2/3 [1 −R(ω)] 1 for a detector volume V with an aspect ratio near unity, and an area A det instrumented with photodetectors. For photon energies ω 1.9 eV or wavelengths λ γ 650 nm, silver has 1 −R(ω) 10 −2 . For higherenergy light, aluminum is better with 1 −R(ω) 10 −1 for all λ > 100 nm. Dielectric coatings can achieve even much lower absorbance values 1 −R 10 −2 in narrow energy ranges over all wavelengths of interest (100 nm-10 µm). High-reflectance coatings can thus allow for a small photodetector-instrumented area A det , of linear size 1/ 1 −R(ω) smaller than the bulk detector size, while keeping η refl ∼ 1.
Transmission efficiency-A Stack detector configuration does not require reflective coatings, since the radiation is focused; in fact, it will typically need antireflective (AR) coatings on the interfaces between the gas (with index of refraction close to 1) and the material separating the different slabs, such as glass or silicon, that will have a substantially higher index of refraction. Graded-index coatings can achieve 0.1% reflection over a broad range of wavelengths, while even better performance may be expected over narrow wavelength ranges with thin-film interference coatings. Bulk absorption by the slab container material must also be considered. Synthetic quartz glasses can easily achieve absorption depths in excess of 1 cm in the wavelength range 180 nm-2.5 µm, while silicon exhibits this property for wavelengths above 1.1 µm, covering our energy range of interest.
Optical thickness-Signal photons can get reabsorbed and subsequently quenched if the gas sample is optically thick. For a low-density gas of twolevel molecules with an electric-dipole transition moment µ 1,0 = | 1|µ e |0 |, the light intensity falls off exponentially as ∝ e −αl , with l the optical path length and α(ω) the attenuation coefficient for light with an angular frequency ω near the molecular transition energy ω 0 : Here, we employ the collisionally-broadened lineshape, though Doppler broadening can also be important at sufficiently low densities. The inverse α(ω) −1 is the mean absorption depth, which can be quite short for |ω − ω 0 | < γ col , as the second line of eq. 100 shows. For this reason, a Bulk detector will rely on fluorescent decays to intermediate levels which are not thermally occupied. It follows from eq. 100 that the associated photons from those decays will travel macroscopic distances, because levels that could resonantly absorb them have exponentially small occupation probabilities p i ∝ e −Ei/T , while off-resonant absorption is highly suppressed by the Lorentzian line-shape factor, so η thick ≈ 1. Even for thin Stack detectors, the on-resonance absorption depth α(ω 0 ) −1 can be smaller than the proposed integrated thickness, D = 1 mm for Phase I and D = 100 mm for Phase II. However, in the regime of eq. 93 where cooperative effects for DM absorption are in action, the radiated photons themselves will also interact cooperatively with the molecular medium. Dicke's seminal work on superradiance showed that due to a phase-matching effect, a photon plane wave can experience coherent scattering in the direction of the incoming wave [5]. This coherent forward scattering in extended volumes was further developed in refs. [18,19], showing that for large-Fresnel-number samples (such as a slab), photons are re-emitted in a narrow "lobes" primarily in the forward direction. Recent theoretical [20] and experimental [21] work shows that this effect is also obeyed for a single photon absorbed by a sufficiently strong dipole transition. Single-photon superradiance in extended volumes is a complicated subject of intense recent study [22,23], and beyond the scope of this work. Further work is needed to show how far a photon will travel in this regime, and to what extent the directional information is preserved after multiple coherent scattering events. This will be crucial for a Phase II version of the Stack configuration.
Photodetector performance-The sensitivity of our proposed setups is ultimately limited by photodetector parameters. The key specifications not only include the quantum detection efficiency η det , but also the dark count rate (DCR), photosensitive area (A det ), energy range (ω min -ω max ), energy resolution (∆ω), timing jitter (∆t), and operating temperature. Reviews on these aspects of single-photon counting detectors can be found in refs. [24,25].
The Bulk configuration requires a large photosensitive area to keep η refl ≈ 1. In a prototype Phase I, we propose utilizing off-the-shelf photomultiplier tubes (PMT) with photosensitive area A det = (30 cm) 2 . These allow for (near-)room-temperature operation at visible and near-infrared wavelengths with detection efficiency η det ≈ 40 %, and DCR ≈ 1 Hz (when they are cooled a few degrees below freezing) [26][27][28][29]. Timing jitter ∆t is sub-nanosecond, aiding in timing-based rejection of environmental backgrounds, while the intrinsic energy resolution is quite poor.
For the Phase II Bulk setup, we assume state-of-theart photodetector arrays based on microwave kinetic inductance detectors (MKID) [30]. These cryogenic photodetectors can be operated as single-photon counters between ultraviolet wavelengths of 100 nm, all the way to mid-infrared wavelengths greater than 5 µm, while retaining energy resolution of order ∆ω ∼ 0.1 eV as well as good timing resolution ∆t ∼ 10 −6 s, with further improvements on the horizon [31]. These devices have essentially no intrinsic dark counts at energies ω a few times above their energy resolution; any non-signal counts must be due to true environmental photons. Their quantum efficiency is already good even for mid-infrared photons (η det > 0.2), and excellent for wavelengths above 500 nm (η det > 0.5).
Crucially for our purposes, it is the first cryogenic photodetector technology to have been multiplexed into arrays of 2 × 10 4 "pixels" [32,33], already yielding photosensitive areas as large as A det ∼ (1 cm) 2 ; larger arrays are already under development [34]. (Most cryogenic single-photon counters rely on the increase in tempera-ture in a superconducting volume from the impinging photon's energy, so any one such volume must necessarily be microscopic.) Anticipating that MKID technology matures even further, we assume a photosensitive area of A det ∼ (10 cm) 2 for the future BII setup volume of V = (2 m) 3 , where it would occupy a 10 −3 fraction of the container area and thus require tuned dielectric coatings with 1 −R ∼ 10 −3 .
In the Stack configuration, the photons can be focused onto a small area of order 10 −6 the transverse area of the molecular container. As such, there is no need for particularly large arrays even in the larger Phase II version, which would require a minimal photosensitive area A det ∼ O(mm 2 ). MKIDs would perform even better as the smaller phase space of the signal photons allows them to make use of microlenses to steer the photons on their inductive elements [31]. Over such relatively small areas, transition edge sensors (TES) [35] and other cryogenic detectors [24,25,36] may also be employed, with similar and potentially better specifications.

C. Environmental backgrounds
Assuming that DCR can be controlled down to the desired levels specified above, we expect three primary sources of background in our energy range: black body radiation (BBR), natural and cosmogenic radioactivity, and cosmic rays.
Blackbody radiation-The molecules in the detector volume have to be kept in the gaseous phase in order to retain their resonant absorption characteristics. For all but a few diatomic molecules suitable for our setup, this corresponds to temperatures larger than 100 K. The accompanying blackbody radiation (BBR) at a temperature T in thermal equilibrium results in an irreducible background photon rate onto the photosensitive area A det of: in a band ∆ω around any frequency ω. Fake counts from BBR will dominate over dark counts in the photodetector whenever η det Γ BBR DCR, cfr. eq. 96. Given the exponential tail in eq. 101, this turnover point occurs at an energy ω that is not very sensitive to the other experimental specifications. We foresee blackbody radiation to be the dominant background for ω < 1.5 eV in any setup at 300 K (BI), and for ω < 0.5 eV in a 100 K detector (BII, SI, SII). These turnover points are clearly visible as the "kinks" in our sensitivity estimates in e.g. fig. 8.
Natural and cosmogenic radioactivity-Background photons in the energy range of interest can also originate from radioactive decays in the surrounding material or even the molecular gas itself. Care must be taken to passively mitigate and, in the Phase II versions of our setups, to actively veto this background. Gamma rays will be the main culprit due to their large penetration length, which is about one and three orders of magnitude larger than that of beta and alpha particles, respectively, at 1 MeV. Gamma rays in the surrounding material can dislodge surface electrons via the photoelectric effect or directly ionize the molecules in the bulk (among other channels), which can cause molecular excitations in the frequency band of interest, or even directly trigger the photodetector.
A later-stage detector proposal should include a detailed detector simulation of these effects; here we provide parametric estimates showing the feasibility of passive and active mechanisms to keep the rate of radioactive decays Γ RD below the inverse shot time or the DCR, whichever is larger. When Γ RD DCR + t −1 shot , radioactive backgrounds are subdominant. The attenuation coefficient for gammas at 1 MeV is α γ ≈ 6 × 10 −2 cm 2 g −1 , giving an attenuation length of about 6 cm in quartz and 1.5 cm in lead for example. The main sources of naturally-occurring radioactivity are 238 U, 232 Th, and 40 K, each giving roughly similar contributions, so we will only estimate contributions from the first. The natural mass-fraction abundance of 238 U is f m ≈ 10 −6 g/g, while its half-life is t 1/2 ≈ 4.5 × 10 9 y. Supposing that the container is a cubic volume with boundary area A δV = 6(2 m) 3 surrounded by a material of quartz's density, we find a typical radioactive decay rate of O(10 4 Hz) in a boundary layer with thickness of order the penetration length. The same estimate for a high-purity lead shield with a specific density ρ Pb ≈ 11g/cm 3 and a f m ∼ 10 −12 mass fraction of 238 U, similar to those used by other experiments [37], gives a rate of Besides naturally-occurring radioactive isotopes with lifetimes on the order of billions of years, there will generally also be trace concentrations of cosmogenicallyactivated radiaoactive isotopes in the molecular medium or the container materials. Even though these contaminants occur in much smaller mass fractions, their decay rates are larger, so that they may compete with radioactivity from 238 U, 232 Th, and 40 K. A particularly dangerous isotope is 14 C, produced by cosmic-ray collisions high up in the atmosphere, and consequently present at a fractional number density of 10 −12 relative to that of 12 C in all organic material not buried for longer than its half-life of t 1/2 ≈ 5730 y. For example, if the molecular medium contains carbon atoms (e.g. CO), then one can expect a cosmogenic radioactivity contribution of: Carbon-containing molecules derived from fossil fuels can have greatly depleted 14 C content. Certain petroleum reservoirs have been shown to contain 10 −18 fractional concentration of the radioactive isotope [38], making this background completely subdominant even at large pressures. We foresee that similar provisions can be taken for other cosmogenically-activated radioactive isotopes. Even a background radioactive rate of Γ RD ∼ 10 −2 Hz, potentially achievable with a high-purity shield and keeping radioactive contaminants at a minimum, poses a huge challenge to a typical DM absorption detector in the ω ∼ eV range. Other detector proposals [39][40][41][42][43][44][45] employ bulk target volumes based on nonresonant absorption onto liquids or solids, and aim to be signal-countlimited at kg-year exposures, requiring powerful active veto methods for large volumes. Our Bulk and Stack configurations each have characteristic properties that provide additional passive mitigation mechanisms not available to the aforementioned nonresonant absorption targets.
The Bulk configuration in scanning mode has three distinct advantages regarding passive radioactive background mitigation. Firstly, by dividing the full integration time over an O(1) bandwidth into many different shots with independent, narrow-band DM response functions, the radioactive background becomes already negligible when Γ RD t −1 shot -a much looser criterion than Γ RD t −1 int . Secondly, one can look for DM signals Γ det smaller than Γ RD even when Γ RD > t −1 shot (as we have implicitly assumed in eq. 96) because there is a natural way to modulate the DM signal rate while keeping the background rate constant. For example, when in neighboring shots one has average total background counts of Γ bckg t shot 1, one is sensitive to signal rates as low as Γ det ∼ Γ bckg /t shot (assuming η det = 1 for simplicity). Such averaging effects and differential response tests are not typically available in most bulk detectors, which do not have a natural way to "turn off" the signal while keeping the background rate constant.
Thirdly, a Bulk detector in scanning mode can be operated at such low density that a typical fast electron only has a small probability P e ex of exciting a molecular state in the frequency band of interest, since that process does not receive a resonant enhancement factor: approximately valid for nσ e ex L 1, and L the typical linear size of the detector volume. Above, σ e ex should be the (velocity-averaged) cross-section for a typical electron to excite molecular states that give rise to fluorescence photons within a detector's bandwidth around the frequency of interest, for which we have chosen a plausible value for an E1-allowed transition in the numerical estimate [46, §148]. Forbidden transitions have even lower σ e ex that scale like α 3 a 2 0 or even higher powers of the fine structure constant. When N e P e ex 1, where N e is the typical number of electrons produced in a radioactive decay, the radioactive background in the band of interest is further suppressed by this factor.
The Stack configuration has the remarkable property that 84% of the signal is emitted in a cone of opening angle 2v 0 ≈ 5.4 arcmin and solid angle ∆Ω = πv 2 0 , a 1.5 × 10 −7 fraction of the full solid angle ∆Ω = 4π. Radioactive backgrounds from the molecular volume can thus likely be ignored entirely, reducing the problem of radioactive contamination only to the photodetector material and mount, a dramatically smaller volume V ∼ O(mm 3 ).
Furthermore, an active veto system consisting of scintillating material and PMTs surrounding the detector volume may be employed. Radioactive decays give rise to many high-energy particles. When they trigger the scintillating material, any "signal" in the photodetector in a short time span around the trigger time in the photodetector can be vetoed. We note that a typical radioactive decay chain releases dozens of primary gammas and betas, as well as a bunch of secondaries, so the trigger efficiency on any one photon or electron need not be high. As long as the photodetector has a sufficiently short jitter time ∆t and the molecular medium a sufficiently fast relaxation rate to thermal equilibrium (through radiative and nonradiative channels), the resulting dead time can be made negligible. Roughly, one requires to leave the duty cycle of the experiment essentially unaffected. The timing jitter requirements are easily satisfied by many orders of magnitude. Figure 6 shows that the detector relaxation is sufficiently fast even for low-lying vibrational levels, if they have dipole-allowed radiative decays. Only for vibrational states that have exclusively E1-forbidden radiative decay channels, does eq. 105 become hard to satisfy. Given the wealth of passive mechanisms on offer to keep the radioactive background under control, the setups under consideration do not require extremely efficient active veto systems, especially in comparison to other proposed experiments in this energy range. The Phase I prototypes can likely forego a large-scale active veto system altogether.
Cosmic rays-Cosmic muons can also make up a considerable fraction of the background. At sea level, the vertical muon flux density is about 70 m −2 s −1 sr −1 [47,48]. The cosmic muon flux is much lower underground, falling to an integrated vertical flux density of 10 −3 m −2 s −1 about 1 km deep in standard rock (2.65 km water-equivalent), and 10 −5 m −2 s −1 at 2 km depth [49][50][51][52][53]. Underground operation in a mine of moderate depth automatically ensures that the cosmic muon background is subdominant to that of radioactivity. Phase I prototypes can likely get away with surface-level operation when outfitted with a modest muon veto system that has a rejection power of 10 2 , one that can possibly work in conjunction with a radioactivity veto.

D. Signal discrimination strategies
From the discussion so far, a signal discrimination strategy emerges. Radioactive and cosmic background events can be identified by the fact that they result in multiple photon counts and ionized electrons. All environmental and detector backgrounds will be distributed over a broad energy range. In contrast, the nearmonochromatic DM signal can only excite transitions at one particular transition energy, and cause a single fluorescence photon at one frequency (or at most a handful of frequencies, if there are decay channels to several intermediate states). Below, we outline that more detailed follow-up studies can unequivocally verify the DM origin of any potential signal, and furthermore determine properties such as mass, spin, interaction type, and 3D velocity with pinpoint precision. The spectacular discrimination power of this type of detector relies on three different "handles": energy response, selection rules, and spatial coherence.
Energy response-The highly resonant response of the detector at low number densities can be used to confine the DM signal to better than 10 −6 fractional frequency precision, and even to perform precision studies of its lineshape. A cryogenic photodetector with energy resolution of ∆ω 0.1 eV can resolve both electronic splittings and vibrational fine structure in most molecules. Temperature modulation will change the relative populations of rotational levels, and thus the absorption rates in the rotational fine structure. By varying the pressure, one can determine how far off resonance the signal is. Scanning with external electromagnetic fields and use of different molecular species-and isotopes of the same species-can be employed to further hone in on the true energy of the signal, i.e. the DM mass.
What is the ultimate energy resolution? Collisional broadening rates can be made arbitrarily small by lowering the pressure, while the fractional radiative linewidth is (γ 0 + i γ i )/ω 0 α 3 ≈ 4 × 10 −7 and typically several orders of magnitude smaller even for dipole-allowed decays. For photon absorption, the Doppler width due to molecular motion often determines the minimal fractional linewidth of the signal at low number densities.
Doppler broadening of the absorption line for a nonrelativistic particle is less pronounced than for a (necessarily relativistic) photon. Using momentum and energy conservation for an inelastic collision with initial (final) molecular velocity v 1 (v 2 ) during which a dark matter particle with mass m, energy ω, and velocity v is absorbed onto an internal molecular state with transition energy ω 0 , we find: The third term on the RHS is quantitatively subdominant, and also does not depend on the molecule's initial velocity. Approximating the dark matter mass with its total energy ω, and expanding in small velocities, we then finally arrive at the condition for absorption: Because the molecules in the gas move at different velocities, they can absorb dark matter particles of differing energies ω. Defining the molecular speed in the dark matter's velocity direction as v 1, ≡ v 1 ·v/|v|, we can express the probability of finding a molecule between v 1, and v 1, + dv 1, as proportional to the Boltzmann-weighted exp{−M mol v 2 1, /2T }dv 1, . Translating this to a normalized molecular transition energy distribution via eq. 107, we find with a fractional width of ∆(v)/ω 0 ≈ 4 × 10 −9 (v/10 −3 ) for H 2 gas at room temperature, and even lower for heavier molecules and/or colder temperatures. This width ∆(v) is a factor of v ∼ 10 −3 smaller than the equivalent Doppler width for photon absorption. In case Doppler broadening does dominate over both radiative and collisional broadening, when ∆(v) γ, then the result of eq. 10 should be convoluted with the line shape of eq. 108 to give: at least for |ω − ω 0 |/∆ not too large. The off-resonance tails are more accurately described by the Voigt profile, a convolution of both the Lorentzian and Gaussian line shapes. The narrow fractional width allows for extremely narrow spectroscopic studies of the signal line shape at low pressures (when collisional broadening can be ignored), and would be a great signal discriminant should anomalous fluorescence be seen in the experiment. A single velocity component of the dark matter field ensemble has an energy ω = m(1 + v 2 /2) in terms of its square velocity in the lab frame. The DM's 3D velocity distribution f (v) can be expected to closely resemble the virialized distribution of eq. 14, yielding a fractional underlying signal frequency width of order v 2 0 /2 ≈ 3 × 10 −7 . Because ∆(v)/ω 0 v 2 0 1, the molecular resonance can resolve the kinetic energy distribution of DM! In this "resolved" regime, the absorption rate as a function of DM mass m is : By tuning ω 0 and at large enough statistics, one could determine the DM mass m, the velocity dispersion v 0 , and the magnitude of the relative velocity |v lab | (using diurnal/annual modulation). One could possibly even discern finer details of the velocity distribution f (v), such as the Galactic escape velocity (not taken into account by the virialized distribution in the second line of eq. 111). Note that the integrand in the second line of eq. 111 is technically only valid for velocities where ∆(v) γ; one should use the full Voigt absorption line shape for the low velocities v where this is no longer a good approximation.
Selection rules-If a near-monochromatic signal were to be detected at some frequency ω near a transition frequency ω 0 , then one can determine the properties (such as the quantum numbers) of the initial and final states. This is possible because small polyatomic systems are simple enough to have been well-characterized both experimentally and theoretically, as is evident from secs. II C and II D.
Any nonthermal SM background must come from interactions with external photons or charged particles, which to leading order interact with a Hamiltonian proportional to the electric dipole moment µ e , giving Ω ∝ |µ 1,0 |. This fact leads to the well-known dipole selection rules for the leading-strength transitions.
DM absorption does not necessarily obey these rules. A scalar DM particle may leave the angular state of the molecule unaffected via a "monopole" transition (∆J = 0), while a pseudoscalar DM particle may cause "spindipole" transitions from a spin-singlet ground state to a spin-triplet excited state (∆S = 0), both at leading order. These processes are normally highly forbidden in small polyatomic systems. Even when the DM particle primarily causes ordinary dipole transitions, one could test if rather than coupling to the dipole moment of electric charge, it instead couples to the dipole moment of e.g. baryon number, lepton number, or any linear combination thereof. Although we do not discuss it in this work, dark matter candidates with spin ≥ 2 would dominantly cause higher multipole transitions.
The modularity of the proposed detector allows for targeted studies with different types of molecules, each having a transition energy ω 0 near a potential candidate signal at ω. Should a DM signal be seen in a broad frequency search, these targeted studies can in principle determine the form of the interaction Hamiltonian δH.
Spatial coherence-The nonrelativistic velocity distribution f (v) leads to a characteristic spatial (and temporal) coherence of the perturbing wave that is imprinted onto the molecular emission via the correlation function in eq. 19, in turn leading to the dramatic focusing effect for a slab-like container discussed in sec. II B. This emission pattern cannot be mimicked by any standard SM background. A robust prediction for such a detector shape is that the center of the emission cones, with offset from the normal direction proportional to v lab (see eqs. 24 and 25), would precess on a diurnal and annual basis due to the Earth's spin and orbit around the Sun, respectively, with known phases, directions, and amplitudes. In addition, these measurements of v lab would have to agree with the magnitude |v lab | derived from eq. 111; likewise, v 0 as determined by the opening angle of the emission cone must agree with the fractional frequency linewidth in eq. 111. These observations can also be compared against astrophysical inferences of the DM's velocity distribution. Finally, a precision line study might reveal additional information not otherwise attainable, such as the existence of DM streams and the relative rotation of the DM halo and the Galactic disk.

IV. DARK MATTER SENSITIVITY
In this section, we present sensitivity projections of the proposed setups to the parameter space of specific DM models, after briefly reviewing the chief interactions of each DM candidate. We have classified the models in terms of spin and parity of the DM boson, starting off with spin-1 vectors in sec. IV A, and then continuing with parity-even, spin-0 scalars in sec. IV B and parity-odd, spin-0 pseudoscalars in sec. IV C. We leave a treatment of DM particles with spin 2 and higher to future work. We summarize the dark matter candidates and couplings, as well as the corresponding transitions they can mediate, in table III.

A. Vectors
We will start with two vector dark matter candidates, namely a hidden photon kinetically mixed with the usual electromagnetic field, and a new photon that couples to baryon-minus-lepton number B − L, both having a Stückelberg mass m γ . We focus on these two cases because of pedagogy (simple comparisons can be made to the interactions of the normal photon), and because they embody simple extensions of the Standard Model that are theoretically consistent up to very high energy scales. Other types of vectors are certainly possible, but most of the models include other states, like Higgs-like scalars or anomalons, whose interactions are often independently constrained; we ignore those theories only for the sake of brevity.
Weakly-coupled, massive vectors can be produced in the early Universe. A natural and calculable relic abundance can arise from vector fluctuations during the inflationary era (if the Stückelberg mass is "on" then), with a present-day relic energy density of where H I is the Hubble scale during the last few e-folds of inflation [54]. We see that the vector can make up all of the DM in the mass range of interest if the Hubble scale is between 10 12 GeV and 10 13 GeV. The field misalignment mechanism, on the other hand, is not effective unless large interactions with curvature invariants are present [55].
Kinetically mixed photon.-If a light vector particle is a low-energy remnant of a sector that is coupled to the SM at high energies, the associated vector field A µ can and will generically have effective operators coupling it to the SM even if none of the SM fields are charged under the new U(1) gauge symmetry [56,57]. The lowestdimensional such operator, one that can be expected to capture the dominant effective interactions with the SM at low energies, is the kinetic mixing term F µν F µν wherein the "hidden" field strength F µν ≡ ∂ µ A ν − ∂ ν A µ couples to the equivalent quantity F µν of the SM electromagnetic field A µ . The strength of this mixing is conventionally quantified by a dimensionless parameter in the Lagrangian: with m γ the hidden photon's Stückelberg mass, and J µ EM = ψ q ψψ γ µ ψ the electromagnetic vector current. The above Lagrangian, written in the "gauge basis", can be made to have diagonal kinetic and mass terms with a field redefinition, transforming it into the so-called physical basis: now rewritten in terms of redefined fields A µ and A µ . Diagonalizing the terms that govern the propagation in vacuum comes at the cost of introducing an interaction of the massive photon state with the EM vector current, suppressed by relative to that of the massless photon. In the physical basis, electromagnetic charges interact with a specific linear combination of fields A eff µ = A µ + A µ , so matrix elements for interactions of A µ with the molecule can be found simply by rescaling those of the photon A µ by a factor of . (This is strictly only true for electric dipole (E1) transitions. For magnetic dipole (M1), electric quadrupole (E2), and even higher-order transitions, whose matrix elements all involve factors of k · x, the matrix elements of a non-relativistic A µ are further suppressed by factors of velocity v = |k|/m γ relative to those of the necessarily relativistic photon A µ .) The propagation of the two vector fields is qualitatively different, however. In vacuum, the A µ = (φ , A ) mass eigenstate obeys the massive wave equation (∂ 2 t − ∇ 2 + m 2 γ )A µ = 0, which means the field can support nonrelativistic solutions as well as a longitudinal modes, which have a vector potential aligned with the propagation velocity v (i.e. A · v = 0). The dark matter state is expected to be well-described by a nonrelativistic, classical solution of this wave equation; a single momentum component of the whole DM ensemble has a vector potential TABLE III. Dark matter candidates classified by their spin, parity, and interaction Hamiltonian δH, along with the types and strengths of transitions they can induce. Transition types considered include those of the electronic (el), vibrational (vib), and rotational (rot) kind; bold face is used whenever the transition type can be E1-allowed and thus used in the Stack configuration. The fifth column contains the diatomic-molecule selection rules on molecular-axis angular momentum projection Λ and inversion i for electronic transitions, and vibrational quantum number v and total angular momentum J for vibrational and rotational transitions. Unless otherwise noted, the total nuclear spin SN and electronic spin Se do not change in absence of spin-orbit coupling. Rabi frequencies Ω are quoted at ω0 = 5 meV, 0.4 eV, 0.8 eV, 8 eV for rotational, ∆v = 1, ∆v = 2 vibrational, and electronic transitions, respectively, using experimentally allowed benchmark values of the coupling at m = ω0. Numerical estimates use FC factors of 10 −2 , ∆J = 1 rotational matrix elements of 1/ √ 3, Re = a0, δe = δe,1 = δe,2 = 1, ∆qB−L = 1, and matrix element estimates from eqs. 47,48,49,53,58,67,74,80. of the form: The full state of the field should be regarded as a mixed state composed out of many of these different velocity components drawn from a probability distribution f (v), each with random phases α v and possibly also random directionsn. We assume the velocity distribution to be close to that of eq. 14, and the directionn of the vector potential to have a coherence time at least as long as 1/mv 2 0 but shorter than the time scale of the experiment. For reasons outlined in appendix A, all of our main re-sults remain valid even when the classical approximation breaks down, at hidden photon masses m γ 15 eV. In this regime, the local dark-matter field occupation numbers become so low that e.g. |A | 2 |A | 2 .
In Lorenz gauge (∂ µ A µ = 0), the hidden electric scalar potential φ can be determined from the relation ∂ t φ = −∇ · A , from which it follows that it is suppressed relative to the hidden vector potential, as typically φ ∼ v|A |. The stress-energy tensor contribution T µν = F µ λ F λν +(1/4)g µν F λσ F λσ +(1/2)m 2 γ A µ A ν from the massive photon state can be evaluated on the solution of eq. 115, and contributes as an effective fluid with nearly zero average pressure, and energy density: Expressed as fields rather than potentials, we see that the hidden magnetic field B = ∇ × A is velocity suppressed relative to the hidden electric field E = −∇φ − ∂ t A , which oscillates with an amplitude of along the direction of the A potential, up to velocitysuppressed corrections. Because of the enormous size of this hidden electric field, which in some ways is equivalent to shining a kilowatt-class "hidden" laser with beam waist of order the detector size of 30 cm (and approaching megawatt power levels for a 3-meter detector size), it is possible to have appreciable event rates in absorptive media even for tiny values of the mixing parameter , as we showed already in fig. 6. A particle ψ with mass m ψ , electromagnetic charge q ψ coupled to the two photons as in eq. 114 has nonrelativistic dynamics dictated by the Hamiltonian where we have ignored spin-orbit coupling, and have assumed in the second line that the particle moves in an electrostatic potential well V = eq ψ φ of other nearby particles, as well as in a massive hidden photon wave of the form 115, with the ellipsis representing terms of O( 2 ) and O(v). An ambient A wave will also interact with any surrounding conductive elements, including the reflective coating of the vapor cell. The resulting screening currents will set the interacting linear combination of vector potentials, A eff = A + A , to near-zero within a hidden photon's Compton wavelength m −1 γ away from the container wall [58]. However, A and A propagate at different speeds, so for a sufficiently large container volume V m −3 γ they will oscillate in and out phase in most of the bulk interior such that we have |A eff | A 0 to a very good approximation.
In a molecule, we can thus separate the full Hamiltonian into the Hamiltonian H 0 from eq. 31 and the interaction Hamiltonian δH = ψ (iq ψ e /m ψ )A · ∇ ψ , with the sum running over all particles in the molecule. Using the identity ∇ ψ = −m ψ [H 0 , r ψ ], we can rewrite the off-diagonal matrix elements-the only ones that matter for transitions-of δH between an initial state |i and a final state |f as: In the second equality, we have made the approximation iω 0 A E , which holds insofar as the hidden photon is on resonance with the transition, i.e. m γ ω 0 . We also assumed that to leading order we can take A(x) and E(x) to be spatially constant, i.e. independent of x, over the spatial extent of the molecule. This approximation is especially precise for a nonrelativistic wave.
To leading order, we can thus think of transitions being caused by the familiar operator of an electric dipole µ e in an effective electric field E eff = E . It can induce transitions in many systems, including vibrational and electronic transitions in diatomic molecules. For vibrational transitions, we can integrate out the electronic motion and regard a neutral diatomic molecule as a spring with charges ±δ e (R) attached to the ends, where we take the charges δ e to depend on the internuclear separation R. (The dependence on R is easy to see: e.g. as R → ∞, it must be that δ e (R) → 0 if the molecule dissociates into neutral atoms.) This function δ e (R) can be calculated from first principles or indirectly measured in absorption spectra. In this simplistic view of the molecule, the transition operator in eq. 120 reduces to ≡δe,1 + . . . .
The first term within curly brackets acts trivially on the vibrational state (it only induces rotational transitions), but the second term can induce ∆v = ±1 transitions as discussed in sec. II D 2, as long as the factor in square brackets δ e,1 is nonzero. In general, δ e,1 is roughly of order the electric dipole moment of the molecule divided by R e for heteronuclear diatomics, and zero by symmetry for homonuclear diatomics. The third term and higherorder terms in the Taylor expansion contribute mostly to higher-harmonic transition matrix elements, and can usually be ignored to leading order, which we shall do so for this discussion. Finally, we find the angle-averaged squared Rabi frequency from hidden-photon dark matter to be for vibrational transitions from the ground vibrational state |v i = 0 to the first excited state |v f = 1 . To get to the second line, we used eq. 117 with ρ γ = ρ DM and eq. 47 with M the reduced mass of the diatomic and ω e the vibrational splitting; the angle-averaged rotational matrix elements are given below eq. 61. For ∆v = 2 transitions, the squared Rabi frequency is reduced relative to that for ∆v = 1 by the factor ω e x e /8ω e , cfr. eq. 50. For electronic transitions, we can repeat the same exercise to find: where | v f |v i | 2 is the Franck-Condon factor from eq. 44, and where we have parametrized the electronic transition moment | χ el f | n er e,n |χ el i | ≡ δ e,2 eR e in terms of the dimensionless number δ e,2 .
In figure 8, we plot the estimated sensitivity for the Bulk I & II configurations in thick blue bands, and for the Stack I & II configurations in thin red bands. For these estimates, we used the Rabi frequency reach of eq. 98 using the configuration parameters of table II in conjunction with the matrix elements of eqs. 122 and 123. We also assumed δ e,1 ∼ 1 and the parametric estimates of eq. 53 for vibrational transitions, and a electronic transition moment with δ e,2 R e ∼ a 0 and a typical FC factor of | v f |v i | 2 ∼ 10 −2 for electronic transitions. Squared rotational matrix elements were conservatively taken to be 1/6.
The BI prototype, the least aggressive design, is seen to be already capable of probing new parameter space for electronic transitions above 1.2 eV. The BII configuration would be a drastic step up in reach for the same electronic transitions, and would also extend the reach to lower masses via ∆v = 2 transitions. Transitions with ∆v = 1 would yield larger absorption rates, but would produce trapped and subsequently quenched fluorescence photons due to the high optical thickness of such a Bulk detector, as we discussed around eq. 100.
The Stack configurations are also optimally operated with molecules exhibiting strongly-allowed electronic transitions at high energies, and ∆v = 2 vibrational transitions between 0.6 eV and 1.2 eV. They can avoid the quenching issue for ∆v = 1 transitions because of their optically thin planar design, and thus have access to larger matrix elements at energies below 0.6 eV. They also pick up less BBR due to their smaller photosensitive area, with the result that even the SI prototype will likely outperform the much larger BII detector at low masses.
In figure 8, we also show a few exemplary sensitivity curves for the proposed experiments. In blue, we show the reach around 2 eV-3 eV from the famous visible-light absorption band B in 127 I 2 , which is both so wide and dense in frequency space that it nearly covers an octave contiguously at standard atmospheric conditions (without need for scanning), making it an attractive molecular candidate for a proof-of-principle prototype experiment. The blue I 2 curve assumes a vapor pressure of 0.25 bar, which occurs in equilibrium around room temperature, and an integration time of 10 5 s, and otherwise Bulk Phase I parameters. Molecular information on iodine was taken from refs. [59][60][61][62][63][64]. Also in blue, at around 0.8 eV, we plot the estimated sensitivity for ∆v = 2 absorption in 1 H 35 Cl at P = 0.25 bar after a single shot of 10 3 s in the Bulk Phase II experiment. We also depict the reach with 12 C 16 O in a Stack Phase I experiment at P = 5 bar, displaying both its infrared ∆v = 1 vibrational transition (see also fig. 7) and its first allowed electronic transition X → A in the ultraviolet. Molecular data on this electronic transition can be found in refs. [65][66][67][68][69][70]. Finally, we depict in red the sensitivity to absorption of hidden photons heavier than 11 eV onto the first three E1-allowed electronic transitions X → A,B,B' in 1 H 2 [71] in the Phase II version of the Stack configuration.
The gray exclusion regions in fig. 8 depict 95% CL constraints on from null observations by Xenon10 [72] of hidden photons from the Galactic DM halo [73] ("Xenon10") and of hidden photons emitted by the Sun [74] ("solar emission"). Although beyond the scope of this work, it would be interesting to work out the detection prospects of this solar emission component in our proposed setups, even though the resonant detector response is likely not optimal for a thermal emission spectrum. Finally, the black line labeled "N eff " indicates an upper bound on derived from the effects that evaporation of hidden photon DM into the photon bath would otherwise have on the effective number of abundant neutrino species in the early Universe [55].
The small number parametrizing the strength of the DM vector interaction is the gauge coupling g. For simplicity of presentation, we assume the new vector does not kinetically mix with the SM photon, and write its Lagrangian as: where the B − L vector current is defined as J µ B−L ≡ ψ q B−L,ψψ γ µ ψ. Analogously to the analysis done for the (kinetically mixed) photon in eqs. 118-120, we can find that the relevant transition operator is that of a B−L dipole moment µ B−L ≡ g ψ q B−L,ψ r ψ coupled to the electric component of the B − L field strength E . We can write the B − L current J µ B−L as the sum of the electromagnetic current J µ EM of eq. 113 and the neutron vector current, as per the charge assignments of eq. 124. Since nuclear motion can be taken to be "frozen" for electronic transitions to leading order in the Born-Oppenheimer approximation, a B − L vector causes the same transition phenomenology as a kinetically mixed hidden photon provided we make the replacement The separability of the B − L current means that we can decompose the dipole moment as µ B−L = (g/e)µ e + gµ n with µ n the neutron number dipole moment. Hence, vibrational transitions in a diatomic molecule are caused by the effective operator: In the first line, we used the same δ e (R) as in eq. 121, and took Z i , A i , and R i to be the atomic number, mass number, and position vector for the ith nucleus, such that the second two terms represent the interaction with the neutron number current. In the second line, we isolated the component of these terms that depends on the internuclear separation R ≡ R 2 −R 1 , and neglected terms acting on the molecular center-of-mass position. So again, we find that, to leading order, vibrational absorption rates of a B − L vector are exactly analogous to those of a kinetically mixed vector, provided we make the replace-ment: This means that one can expect any E1 transition to be sensitive to absorption of a B − L vector, barring accidental cancellations. In addition, even electric-dipoleforbidden transitions, with δ e,1 = 0, may give appreciable absorption rates. For example, the molecule 1 H 2 H has no electric dipole transition moment (by symmetry), but has In figure 9, we plot the reach estimates for our four proposed configurations. They are exactly analogous to those of a kinetically mixed photon with the replacements of eqs. 127 and 129, except now a Bulk detector can potentially also look for E1-forbidden ∆v = 1 transitions with low quenching rates. The exclusion regions from the hidden photon apply in exactly the same way, modulo the rescaling of eq. 127. In addition, a B − L vector mediates a finite-range "Coulomb" force between electrically neutral bodies: the electrons and proton charges cancel each other out, leaving a like-charge repulsive force between neutron pairs. Null results from short-distance gravity tests [75] place significant constraints at low masses, indicated by the gray exclusion region labeled "fifth forces" in fig. 9.  fig. 8, with the addition of an exclusion region coming from searches for short-range fifth forces. We also indicated in dashed green the B − L gauge coupling that would yield a Coulomb force between nucleons with a strength equal to that of gravity.

B. Scalars
In this section, we first summarize the basics of scalar dark matter, which is analogous to that of vector dark matter in sec. IV A and nearly identical to that of pseudoscalar dark matter in sec. IV C. Afterwards, we will demonstrate how the proposed experiment can be sensitive to scalar couplings to electrons, photons, quarks, and gluons If a light scalar (or pseudoscalar) field exists in the spectrum of the theory, then cold ensembles of scalar particles are naturally produced in the early Universe via the field misalignment mechanism [76][77][78]. As in the massive vector case, an abundance of nonrelativistic, weakly coupled scalars form an effectively pressureless fluid, and are thus excellent dark matter candidates (see [79] for a review of the relevant cosmology). With a Lagrangian of the form and sufficiently weak couplings of φ in δL, scalar dark matter is expected to be a mixed-state superposition of plane waves of the form drawn from a velocity distribution f (v) like that of eq. 14 with random phases α v . The amplitude φ 0 can be thought of as the typical magnitude of the field oscillation; more precisely, we take it to be φ 0 = 2 φ(t, x) 2 v,α . The field configuration then carries an energy density of approximately ρ = m 2 φ φ 2 0 /2. If it makes up all of the local DM energy density, the field amplitude is expected to be: where we have normalized the amplitude in Planck units (G N is Newton's gravitational constant). We will use the same notation for the dimensionless fieldφ ≡ √ 4πG N φ. A light, parity-even scalar may couple to parity-even matter operators (in which case the scalar is often called a modulus or dilaton-like field) of the form: Above, m e is the electron mass, ψ e is the electron field, F µν is the electromagnetic field strength, γ q is the anomalous dimension of the quark field ψ q (up ψ u , down ψ d , and strange ψ s ) with mass m q , β 3 and g 3 are the QCD beta function and gauge coupling, and G A µν is the QCD field strength. The dimensionless couplings d me , d e , d mi , d g parametrize the strength of the leading linear couplings of φ to electrons, photons, quarks, and gluons, respectively, here written in an effective Lagrangian at a scale just above the QCD confinement scale Λ 3 . In most cases, one can ignore the effects of other higher-dimensional operators of SM fields, and couplings quadratic and higherorder in φ. (One exception is when the linear couplings are absent or highly suppressed, e.g. via a parity symmetry under φ → −φ. In this case, all of our results apply with a straightforward replacement of φ ↔ φ 2 , ω 2m φ , and an appropriate field rescaling, while constraints from other experiments typically change qualitatively.) We have written the couplings of the scalar field as a low-energy effective theory at the GeV scale, only having included the most relevant operators while remaining agnostic about the theoretical origins of the scalar in the ultraviolet. The simplest UV completion of the couplings in eq. 133 is that of the linear Higgs portal [80] with b a dimension-1 coupling and H the SM Higgs field. At energies far below the Higgs mass, φ inherits part of the couplings of the Higgs to other SM fields due to the small mixing term that eq. 134 induces. Fermion couplings, including to electrons and quarks, are of order: with m h ≈ 125 GeV the Higgs boson mass, while couplings to gauge bosons, d e and d g , are suppressed [80]. We indicate on the plots in figs. 10 and 11 the induced couplings for the Higgs portal model for b < m φ . (Note that the relative couplings to light SM fields in the Higgs portal model are correlated; in particular, their ratios are fixed, e.g. as in eq. 135.) Larger couplings naively would destabilize the scalar potential, barring any other mechanism or fine tuning. Light scalar fields with parity-even couplings can also arise in theories with an extended gravitatational sector, such as the dilaton in string theory [81,82] or a radion in a theory with extra spatial dimensions [83][84][85]. In addition, theories with spontaneously broken supersymmetry and/or flavor symmetries usually abound in light moduli fields [86]. All of the above examples have extra structure beyond that indicated in eq. 133. If new fields associated with this new dynamics come in at a scale Λ, then φ generically receives a mass-squared correction of order: where y e and y q are the electron and quark Yukawa couplings to the Higgs field, and we have assumed Λ m h . Without accidental cancellations, one would expect m 2 φ ∆m 2 φ , an inequality we plot for each of the couplings individually in figs. 10 and 11 for Λ = 10 TeV, an energy scale not yet directly explored in collider physics. We stress that parameter space above this "natural" region is by no means excluded or unattainable. The scalar may be regarded as a composite particle (not unlike the pion in QCD) at a much lower scale in some theories; see [87] for a notable radion construction along these lines. Alternatively, there may be anthropic pressures for tuning the mass of the DM particle.
In a CP-violating background such as a QCD vacuum with nonzero θ angle, the QCD axion a will also pick up parity-even couplings to mesons and nucleons [88], which can be expressed in terms of an equivalent quark coupling bounded to the interval: with f a the QCD axion decay constant. The upper bound comes from null results in searches for a neutron electric dipole moment, constraining |θ| 10 −10 , while the lower "bound" is the minimum natural size predicted in the Standard Model. At higher loop order, one would also expect the other parity-even couplings, but here we focus on the quark coupling only. The effective Lagrangian in eq. 133 is written in such a way [89] that low-energy masses and couplings have the simple functional dependence on the background value of φ: We will often focus on the symmetric combination of the light quark masses, which has a dependencê withm q ≡ (m u + m d )/2 and dm q ≡ (d mu m u + d m d m d )/(m u + m d ). A neutral atom with nucleon number A and atomic number Z has a mass M that scales nearly linearly with Λ 3 but also functionally depends on the pion mass (and thus the sum of quark masses m u + m d ) and the fine-structure constant due to binding energy effects, as well as the electron mass: where we left out subdominant terms coming from e.g. the strange quark mass dependence of the nucluear mass, or its dependence on the light-quark mass difference m d − m u . The "dilaton charges" Qm q , Q e and Q me have been worked out in ref. [89], and roughly obey the following empirical formulae across the periodic table: Spatial and temporal field oscillations in φ such as those in eq. 131 give rise to fractional variations in e.g. the electron mass with amplitude given by the coupling constant d me times the field amplitude in eq. 132. The vibrational Hamiltonian of a diatomic molecule (having integrated out electronic motion) in the presence of a modulus field is: We defined by Qm q ,eff ≡ (Qm q ,1 M 2 + Qm q ,2 M 1 )/(M 1 + M 2 ) the effective dilaton quark charge of the reduced mass M ≡ M 1 M 2 /(M 1 + M 2 ) of the diatomic molecule.
In the second line, we have only keptφ terms with offdiagonal matrix elements, and see that all monopole operators listed in eq. 46 are generated. The first interaction from δH 0 I is the quantum-mechanical operator corresponding to the classical effect of fractional oscillations in the equilibrium size of atoms, upon which the modulus dark matter searches of [90,91] are based. The other two terms are reminiscent of classical parametric resonance, in that they primarily cause ∆v = 2 transitions when the driving field oscillates at twice the harmonic oscillator frequency.
To leading order, we can ignore the spatial dependence in φ(t, x) in the matrix elements of the operators in eq. 147. However, we have neglected a tidal force on the diatomic molecule that comes about from the fieldinduced spatial variation of the nuclear masses: Here, N = 1, 2 runs over the nuclear labels, and j runs over the dilaton charges, the two most important of which are listed in eqs. 144 and 145. This spatial variation leads to a dipole Hamiltonian that acts on the internal state of the molecules: with ∆Q j ≡ Q j,2 − Q j,1 the difference in atomic dilaton charges. The ellipsis indicate terms acting on the centerof-mass degrees of freedom. Heteronuclear diatomics with significantly different dilaton charges may thus experience dipole vibrational transitions as well. Amusingly, the interaction in eq. 152 is the spectroscopic analog of the differential force macroscopic bodies experience in a dark matter background of very light moduli [92,93]. Finally, modulus couplings to electrons and photons can cause monopole electronic transitions of the type anticipated in eq. 55. In the presence of a modulus field, and with the nuclear motion frozen, the electronic part of the Hamiltonian in eq. 31 becomes: Extracting the φ-dependent parts, and using the identity of eq. 58, we find a transition operator whose transition matrix elements have to be estimated or calculated from first principles. As a first approximation, we use the NDA estimate of eq. 58. The spatial variation in the electron mass m e [φ(t, x)] can also lead to dipole transitions via the operator ω < 0.6 eV 0.6 eV < ω <  For δH 1 I , we assumed differential dilaton charges of ∆Qm q ∼ 1/30, ∆Qe ∼ 1/300, and ∆Qm e ∼ 1/4000. For electronic matrix elements, we utilized eq. 58 for δH 0 IV and eq. 67 for δH 1 II . We also assumed v lab ∼ 10 −3 and φ makes up all of DM such that eq. 132 holds. whose transition matrix elements are suppressed by a factor of v lab /α relative to those of the operator in eq. 154 however. It is clear from the form of eq. 155 that we may view this interaction as that of an effective electric dipole operator with [eE] eff = √ 4πG N d me m e ∇φ(t, x). At subleading order in the Born-Oppenheimer approximation, modulus-induced spatial variation in the nuclear mass can also cause electronic transitions and thus extend the sensitivity to d g and dm q to higher masses, but with greatly reduced transition amplitudes.
In figure 10, we show the reach of the proposed setups to the scalar electron coupling d me as a function of the scalar mass m φ , while the three panels of fig. 11 show the same for the scalar couplings to photons, gluons, and light quarks, respectively. At any point in those parameter spaces, there are several operators contributing to absorption, although one is typically dominant. In general, the Bulk and Stack setups achieve highest signal rates for different operators even at the same point in parameter space, because the Stack setup is insensitive to E1forbidden transitions, so it can never probe e.g. monopole operators. We have summarized in table IV which operators of eqs. 148, 149, 150, 152, 154, and 155 drive the estimated optimal sensitivity reach in figs. 10 and 11, and our assumptions for the matrix element sizes. We broke down the parameter space into three regions: ω 0.6 eV (where both ∆v = 1 and ∆v = 2 transitions can occur), 0.6 eV ω 1.2 eV (only ∆v = 2 transitions), ω 1.2 eV (only electronic transitions).
Short-distance tests of the inverse-square law for the gravitational force between neutral macroscopic bodies strongly constrain the low-mass end of the parameter space of interest [75]. Given typical dilaton charges of neutral bodies as listed in eqs. 144 and 145 for the quark and electromagnetic couplings, and assuming Q me ≈ m e /2m p for the typical electron-mass dilaton charge, we recast these "fifth-force" tests as constraints on the individual couplings d j . Regions excluded at 95% CL are shown in gray in figs. 10 fig. 10, except the cooling bound on de comes from cooling processes in the Sun rather than red giants. We also show the experimentally allowed region for the scalar coupling dm q of the QCD axion.
the exact analogue of eq. 132. In other words, the neutron EDM is no longer zero (or constant), but takes on the field-dependent value [100]: with a coefficient d θ ≈ 2.4 × 10 −16 e cm carrying a 40% fractional uncertainty [101]. The QCD axion operator of eq. 157 also gives irreducible, low-energy contributions to the pseudoscalar operators of eq. 156. In particular, one has G aγγ ≈ −1.92(4)α/2πf a primarily from mixing with the pion, as well as G app ≈ 0.47(3)/f a , G ann ≈ 0.02(3)/f a , and a two-loop suppressed G aee [102] for the purely hadronic QCD axion in the KSVZ benchmark model [103,104].
Axions also generically emerge out of string theory with exponentially suppressed masses [107,108], and other light axion-like particles naturally arise as pseudo-Nambu-Goldstone bosons of global symmetries broken at a high scale f a , with low-energy masses and couplings typically scaling inversely proportional to f a .
The first set of derivative interactions of a with the proton p, the neutron n, and the electron e in eq. 156 lead to the nonrelativistic interaction Hamiltonian (for single particles) with f = p, n, e. The first term in square brackets has a trivial action on the molecular wavefunction, and can only generate spin flips; it is the basis for cosmic axion search proposals at much lower masses [100,[109][110][111], as well as searches for axion-mediated monopole-dipole and dipole-dipole forces [112][113][114][115]. The second term in square brackets can excite (spin-)dipole transitions in molecules at much higher energies, and is the one we focus on here. By NDA in nonrelativistic molecules, it can be seen that the pseudoscalar coupling to photons, the second term in eq. 156, causes transitions subleading in strength as compared to the nuclear and electronic coupling in typical models, so we will not discuss it here. The coupling to protons and neutrons generates a coupling to nuclei of the form δH = G aN N (∂ t a)σ N · −i∇ N M N , where we assume for brevity that G aN N = G app ∼ G ann is the same for every nucleus with spin. In general, the coefficient will depend on the nuclear species N as G aN N = c N,p G app + c N,n G ann with c N,p and c N,n coefficients of O(1), but we will ignore this complication. In a diatomic molecule with a spinless nucleus N = 2, and in an unpolarized spin state for the N = 1 nucleus, we find the following square Rabi frequency for ∆v = 1 vibrational transitions: To get to the second line, we used the vibrational amplitude of eqs. 71 and 74, and assumed that the pseudoscalar is all of the DM. The averaged rotational amplitudes are given below eq. 61. Generalizations to amplitudes with multiple nuclear spins, polarized spin states, or nonharmonic ∆v = 2 transitions are straightforward. The irreducible neutron EDM operator of the QCD axion also gives rise to the operator with E(R N ) the internal electric field of the molecule evaluated at the nuclear position R N . The second equality follows from E(R N ) = −∇ N V (R N )/(eZ N ) with V the potential energy terms of H 0 in eq. 31 and using canonical commutation relations. Famously, this operator has vanishing diagonal matrix elements due to Schiff's theorem [116], but its off-diagonal matrix elements do not vanish. Assuming that the nuclear EDM is similar to that of the neutron, d N ∼ d n , and again starting with a diatomic with a single, unpolarized spin σ 1 , we find a square Rabi frequency: We see that the nuclear pseudoscalar operator from eq. 161 generates exactly the same transitions as the EDM operator of eq. 163. The latter's Rabi frequency is always smaller though, by the fraction d θ M 1 /Z 1 e 2d θ m p /e ≈ 2.2 × 10 −2 for G aN N f a = 1. Since its only effect is to give a subleading contribution to transitions already caused by the G app and G ann couplings-which are always generated as well-we will not make separate sensitivity projections for the EDM operator. If the spindipole transition acts trivially on the nuclear spin state, then the transition can be E1-allowed; if the transition is associated with a combined spin flip in one or more nuclei, then it is always strongly forbidden. Electronic transitions may be excited via the pseudoscalar coupling to electrons in a similar way: 1|δH|0 = iG aee (∂ t a)ω 0 1| Z1+Z2 n=1 σ e,n · r e,n |0 ; (165) We have used the matrix element identity of eq. 80 and defined its strength with δ e,2 as in eq. 123. We have notationally suppressed spin degrees of freedom. In a magnetic field, ∆Σ = ±1 spin-flip transitions lines will receive Zeeman splittings relative to the spin-preserving ones ∆Σ = 0, of size given in eq. 90. This Zeeman tuning can be used in the Bulk configuration to achieve efficient and uniform frequency coverage, as all lines receive a uniform shift by the same absolute amount. In absence of significant spin-orbit coupling, ∆Σ = 0 emission rates are suppressed, so this scanning method is less suitable for use in the Stack setup with small molecules.
In figure 12, we plot gross sensitivity projections to the nuclear coupling G aN N with the four proposed setups for the ∆v = 1 vibrational transitions calculated in eq. 162, as well as for the first higher-harmonic ∆v = 2 absorption. Figure 13 contains similar curves for the electronic transitions of eq. 166 induced by G aee . We also plot the typical mass-coupling relation-f −1 a /2 < G aN N < f −1 a with m a and f a related as in eq. 158-for the QCD axion. Barring accidental fine tuning in two separate couplings-G app and G ann -any QCD axion model must lie in or above the green band in fig. 12. Figure 13 depicts by green bands the relation between the electron coupling G aee and the axion mass m a for the DFSZ and KSVZ benchmark models of the QCD axion. The KSVZ band is roughly the smallest electron coupling an untuned model can exhibit. In these sensitivity plots, we make the minimal but optimistic assumption that the axion or axion-like particle makes up all of the dark matter energy density. If the QCD axion is to make up all of the DM, most cosmological histories would predict or strongly favor axion masses below 1 meV if the cosmological abundance arises due to the misalignment mechanism [117], and below 4 meV (with large uncertainties on this bound) if the axions can be produced from decays of topological defects [117][118][119][120][121][122]. In a standard thermal history, the QCD axion abundance scales as Ω a ∝ m 1.19 a in the dilute instanton gas approximation. Therefore, the standard production mechanisms would predict a fractional axion DM abundance of ρ a /ρ DM ∼ O(10 −3 ) if the QCD axion were to exist with a mass m a ∼ 1 eV. One could conceivably construct models wherein these large-mass axions do constitute all of the DM, but we will not attempt to do so here. We note that even if the QCD axion only makes up such a small subcomponent of the DM, the SII setup would still be capable of detecting this component between 0.4 eV and 1.2 eV, as the coupling sensitivity only scales as the square root of the energy density δG aN N ∝ √ ρ a .
The dominant astrophysical bounds on the pseudoscalar couplings to nucleons and electrons come from supernova and white dwarf cooling, respectively. The observed duration of the neutrino burst originating from SN1987a indirectly constrains other efficient cooling mechanisms like the emission of light pseudoscalars. A simple energy loss argument was long thought to set a rough bound of G aN N 2.5 × 10 −9 GeV −1 [123]. A more detailed analysis, correcting the omission of several physical effects in the early literature and folding in progenitor uncertainties, finds a much weaker robust exclusion bound of G aN N 1.7 × 10 −7 GeV −1 [124], which is the one plotted in fig. 12. Observations of drifts in light-pulsation periods in white dwarfs provide an indirect measure of their cooling rate; the most stringent constraint, shown in fig. 13, is G aee < 2.5 × 10 −10 GeV −1 at 95% CL [123,125,126].

V. DISCUSSION
We have presented a DM detection scheme based on resonant absorption in molecules. Our proposed setup is sensitive to a wide variety for bosonic DM candidates with masses between 0.2 eV and 20 eV, including axions, dark photons and moduli, and can achieve several  12. Reach for the pseudoscalar coupling GaNN to nucleons as a function of pseudoscalar mass ma, for the configurations BII (thick blue band), SI, and SII (thin red bands). The green band indicates the typical mass-coupling relation between ma and GaNN ∼ 1/fa for the QCD axion. The gray exclusion region is a "robust" bound set by a combination of the neutrino burst duration of SN1987a, and null observations of axion scattering events in waterČerenkov detectors after the same supernova. orders of magnitude improvement in coupling on current limits over the energy range under consideration. The detector concept may be regarded as a hybrid between low-energy macroscopic oscillators-circuits, cavities, mechanical resonators, electron and nuclear spin resonance systems-and high-energy absorption onto highdensity targets. On the one hand, molecules can be viewed as some of the smallest and highest-frequency electromechanical resonators that exist in Nature. Each molecule is tiny and extremely weakly coupled to DM, but a gas of them contains an enormous number of es-sentially identical copies, boosting the signal rate. On the other hand, our setup is in some ways similar to proposals aiming to look for DM absorption onto bulk targets [39][40][41][42][43][44][45]. The crucial difference is that instead of looking for DM-induced excitations in a continuumsuch as the free-particle continuum, a conduction band, or a phonon spectrum-here we advocate looking for DM absorption into a resolved discretuum of states with a long lifetime T 1 and phase-coherence time T 2 .
A resonant approach comes with numerous advantages, and one major challenge: frequency coverage. Resonant absorption onto a transition line at ω ω 0 only yields appreciable event rates in a narrow bandwidth of order ∆ω ∼ 1/T 2 around any one nominal transition frequency. Fortunately, small polyatomic molecules are multimode resonators due to their wealth of electronic levels, each with their own vibrational and rotational fine structure. In a narrow band around each of the lines in this "forest", the absorption cross-section is resonantly enhanced by a factor of ω 0 T 2 . This enhancement in the absorbing power allows for excellent DM sensitivity at a discrete set of energies with exposures as small as 10 −5 kg year.
Most of the advantages of our proposal boil down to one key feature of molecular spectroscopy: control. Small polyatomics exhibit a resolved discretuum of states with spectra, dynamics, selection rules, optical properties, and response to external variables that are well understood both experimentally and theoretically. As such, a gas of molecules has many "knobs and handles" to control the susceptibility to any DM absorption signal, unlike most bulk absorption targets. Even if a candidate signal is first seen in a different detector, our proposed setup is ideally suited to perform precision follow-up studies, so in this sense it is complementary to other detectors.
Our detector setup allows for excellent background rejection and signal discrimination. Environmental backgrounds can be naturally suppressed due to the low density of the target material, and on otherwise forbidden transitions. The differential energy response of the molecular sample allows for detectable signal rates lower than background rates upon averaging. In addition, background events already become negligible if they occur once per "shot", rather than once over the lifetime of the experiment. Active veto systems may also be employed, as the overall detector has a fast response and relaxation time. Finally, the spatial coherence of the DM particles leads to a dramatic focusing effect of the signal photons in the Stack detector, offering a factor of greater than one million in directional isolation from environmental backgrounds.
Should a signal be found, the combination of great intrinsic energy resolution and energy response control with external variables-pressure, temperature, electromagnetic fields, molecular species and isotope-means that the DM mass can be pinpointed with extreme accuracy and precision, easily resolving even its line shape. By using a variety of molecules with a transition line near this energy, detailed information can be gleaned about the DM's interaction Hamiltonian and selection rules. A dedicated array of Stack detectors could be built to exploit and learn about the kinematic and directional properties of DM. Any signal can be unambiguously identified to have a DM origin. Moreover, a positive signal in this energy range would open up a field of observational DM astronomy given the setup's resolution to both the energy and momentum vector of the DM particles.
We see ample research opportunities for the near and far future. First and foremost, it is imperative that prototypes similar to the proposed Phase I setups get off the ground to demonstrate the feasibility of the experimental strategies outlined in this work. Experimental research and development should include: optimization of MKIDs to deal with isotropic fluorescence photons in the Bulk configuration or other forms of low-noise calorimetric photon detection; manufacturing methods for a physical or artificial set of slabs in the Stack configuration; identification of optimal detector elements, including shield and container materials, and reflective and anti-reflective coatings. On the theoretical front, a significant effort should be devoted to mapping out which molecules are most promising to cover the two decades of energy with their discretuum of transition lines. In addition, optical thickness issues and the effects of multiple forward scatterings of emission photons need to be studied in the context of a Phase II version of the Stack configuration. There could be other uses of our detector concept, such as inelastic scattering of DM, or detection of particles from other cosmic or astrophysical sources. We have considered absorption via a large-but incomplete-set of DM candidates and couplings; notable omissions include the pseudoscalar coupling to photons, and spin-2 DM candidates.
Molecular-gas-based DM searches define an entirely novel class of detectors without the need for novel or exotic material. We have shown that resonant absorption onto molecular transitions is a promising avenue for detection of well-motivated DM candidates, and a new application of molecular spectroscopy and low-noise photodetectors in fundamental physics. In forthcoming work, we will study how variants on the techniques presented here can potentially extend the reach down to DM masses as low as 10 −4 eV.
where we have taken the local DM energy density ρ ≈ 0.4 GeV/cm 3 and velocity dispersion v 0 ≈ 235 km/s. We thus find that N is close to unity for the upper end of the parameter space under consideration in this work. The semiclassical approximation can be expected to be reasonably accurate at low masses, while for m 15 eV the number density becomes so low that a more appropriate representation of the DM interaction is that of individual DM particles impinging on molecules with a small absorption cross-section. Below, we sketch out a fully quantized treatment that encompasses both regimes. The main result will be that all semiclassical results presented in this paper are valid-even when N 1-for the integration times under consideration.
We can write the interaction Hamiltonian a molecular system with a bosonic field mode of energy ω and annihilation (creation) operator a (a † ) as The molecule is approximated by a two-level system consisting of states |0 , |1 with energy splitting ω 0 and interaction-picture annihilation (creation) operators be −iω0t (b † e +iω0t ) as before. In eq. A2, we have (for now) ignored spatial dependence of the local interaction, and absorbed all phases and other numerical constants intõ Ω. Generalizations to interactions with multiple bosonic field modes is straightforward.
We are interested in calculating the molecular transition rate of the process |0 → |1 given an initial darkmatter state |DM of the bosonic field mode under consideration. To this end, we compute the partial rate amplitudes 1; n| t 0 δH (t )dt |0; DM for a combined transition of |0 → |1 in the molecule and |DM → |n in the DM field mode, as an absorption event will in general also affect the state of the dark-matter field. We take |n to be members of an orthonormal basis (e.g. Fock states). The expected absorption probability P abs per molecule is found by summing over the squares of all partial amplitudes with different DM final states |n , which to first order in perturbation theory then gives, at short times: To get to the third line, we have used the fact that the offdiagonal matrix element of b † is unity, and that n |n n| is the unit operator by construction. Matching to the semiclassical result of eq. 2 can be done with the identification: (A4) In the second equality, we defined the number operator N ≡ a † a, and wrote the expectation value as an operator trace weighted by the DM density matrixρ DM (not to be confused with the DM energy density ρ) to allow for the possibility of mixed states.
The identification of the expectation value | Ω| 2 N with a semiclassical perturbation of strength Ω 2 does not quite capture all of the physics, for the operatorN has quantum fluctuations of its own. To see this, one could compute the moments ofN , and see that its variance N 2 − N 2 becomes large compared to its squared expectation value N 2 at low mode occupation numbers. However, these fluctuations-due to the "particle discreteness" of the DM state-average down to negligible levels over the long integration times and macroscopic detector volumes under consideration in this work.