Light Dark Matter in Superfluid Helium:Detection with Multi-excitation Production

Light Dark Matter in Superfluid Helium: Detection with Multi-excitation Production

Simon Knapen Theory Group, Lawrence Berkeley National Laboratory, Berkeley, CA 94709 USA Berkeley Center for Theoretical Physics, University of California, Berkeley, CA 94709 USA    Tongyan Lin Theory Group, Lawrence Berkeley National Laboratory, Berkeley, CA 94709 USA Berkeley Center for Theoretical Physics, University of California, Berkeley, CA 94709 USA    Kathryn M. Zurek Theory Group, Lawrence Berkeley National Laboratory, Berkeley, CA 94709 USA Berkeley Center for Theoretical Physics, University of California, Berkeley, CA 94709 USA
July 15, 2019
Abstract

We examine in depth a recent proposal to utilize superfluid helium for direct detection of sub-MeV mass dark matter. For sub-keV recoil energies, nuclear scattering events in liquid helium primarily deposit energy into long-lived phonon and roton quasiparticle excitations. If the energy thresholds of the detector can be reduced to the meV scale, then dark matter as light as MeV can be reached with ordinary nuclear recoils. If, on the other hand, two or more quasiparticle excitations are directly produced in the dark matter interaction, the kinematics of the scattering allows sensitivity to dark matter as light as keV at the same energy resolution. We present in detail the theoretical framework for describing excitations in superfluid helium, using it to calculate the rate for the leading dark matter scattering interaction, where an off-shell phonon splits into two or more higher-momentum excitations. We validate our analytic results against the measured and simulated dynamic response of superfluid helium. Finally, we apply this formalism to the case of a kinetically mixed hidden photon in the superfluid, both with and without an external electric field to catalyze the processes.

cftsecindent cftsecnumwidth cftsubsecindent cftsubsecnumwidth

I Introduction

Weakly Interacting Massive Particles (WIMPs) with a mass of GeV have been one of the leading dark matter (DM) candidates for the past few decades. However, recent null results in direct detection and collider experiments now provide strong motivation to extend the scope of our models and searches as much as possible. In addition, theoretical advances have shown that there are a variety of models for sub-GeV dark matter that are only now beginning to be explored. Such dark matter may reside in a low mass hidden sector (or “hidden valley”) at the MeV-GeV scale Strassler:2006im (), with either strongly or weakly interacting dynamics Boehm:2003hm (); Pospelov:2007mp (); Hooper:2008im (); Feng:2008ya (). These particles can be invisible to production at colliders, but give rise to large scattering cross-sections in direct detection experiments. They are moreover well-motivated in Asymmetric Dark Matter (e.g. Kaplan:2009ag ()), supersymmetric hidden sectors Morrissey:2009ur (); Baumgart:2009tn (), and SIMP dark matter Hochberg:2014dra (), to name a few.

The theoretical progress in identifying sub-GeV dark matter has been accompanied by effort to experimentally probe such light dark matter Alexander:2016aln (). Among the various ways to detect dark matter, existing direct detection experiments have traditionally focused on nuclear recoils from WIMPs, with typical recoil energies of 10-100 keV. Rapid progress in recent years has thus produced strong limits on DM-nucleon scattering in the 10-100 GeV mass range Aprile:2012nq (); Akerib:2015rjg (); Agnese:2014aze (); Tan:2016zwf (), tightly constraining many well-motivated models of WIMP dark matter. To improve sensitivity to lower mass DM, a number of these experiments have successfully developed techniques that lower the nuclear recoil thresholds below keV. This has been implemented for example in CDMSlite Agnese:2015nto () and CRESST Angloher:2015ewa (), which are sensitive to GeV-scale dark matter.

For a given deposited energy, sensitivity to lighter dark matter can be obtained by scattering from electrons, rather than nuclei. This is because in elastic scattering with the target at rest, the deposited energy is , where is the target mass and the momentum transfer is given by the dark matter velocity and the dark matter-target reduced mass . The first effort in this direction utilized an electron ionization process in XENON10, deriving a constraint on electron interaction cross-sections for DM heavier than 10 MeV Essig:2012yx (). For this mass, the DM possesses the minimum kinetic energy needed to ionize an electron from xenon, 12 eV. In the future, SuperCDMS may have sensitivity to MeV-scale DM, on account of the smaller eV excitation energy set by the band gap of the semiconductor Essig:2011nj (); Lee:2015qva (); Essig:2015cda (). (SuperCDMS may also probe unexplored parameter space for light bosonic DM with eV-keV mass through an absorption process Hochberg:2016sqx (); Bloch:2016sjj ().) Other small gap materials may also make good targets for MeV-GeV mass dark matter in scattering, most notably graphene Hochberg:2016ntt (), giving access to directional information, and crystal scintillators Derenzo:2016fse ().

To reach DM lighter than an MeV, new ideas are needed. The first proposal sensitive to keV-scale DM considered superconductors Hochberg:2015pha (); Hochberg:2015fth () for DM-electron scattering. A conventional superconductor has a small 0.3 meV electron band gap and a large electron Fermi velocity, ; these two facts combined kinematically allow access to keV mass DM (carrying a meV of kinetic energy). It was also shown that these targets have a remarkable sensitivity to bosonic DM in the meV-eV mass range via absorption on electrons, followed by phonon emission Hochberg:2016ajh (). Aside from the meV electron band gap, superconductors have another property which can allow for detection of small energy depositions; a DM scattering that breaks a Cooper pair will give rise to long-lived quasiparticle excitations (which behave very much like an electron). In a very clean superconductor, excitations created in the bulk can then be detected in sensors at the surface of the target. Among the experimental challenges to implementing this idea, it is necessary that the energy resolution be improved significantly, down to the meV scale.

In this paper, we turn to a new proposal to detect keV-MeV scale DM via nuclear recoils in superfluid helium, first discussed in Ref. Schutz:2016tid (). Similar to the superconductor target, the low-energy degrees of freedom in superfluid helium are long-lived quasiparticles. These quasiparticle excitations (called phonons, rotons and maxons) are collective modes in the fluid, analogous to sound waves in the long-wavelength limit. These modes are produced by nuclear scattering, and have been extensively probed by neutron scattering experiments on superfluid helium. Since a large fraction of the deposited energy in a low-energy nuclear scattering is converted to phonons and rotons, it may then be possible to detect dark matter as light as MeV via regular nuclear recoils if experimental thresholds can be lowered to meV. This is because MeV mass dark matter deposits meV of energy in a nuclear recoil process, but this energy is amplified by meV through the evaporation of the excitation at the surface of the superfluid. For a discussion of experimental aspects of a liquid helium detector, and possibilities for detecting the phonons and rotons, see for example Refs. Guo:2013dt (); mckinseyslac ().

The idea of Ref. Schutz:2016tid () was to probe lower mass DM, in the keV-MeV range, by taking advantage of multi-excitation production in superfluid helium. For these low masses, which have correspondingly small momentum keV, the DM couples directly to the collective quasiparticle modes. However, the kinematics prohibit the creation of a single excitation with energy above a meV. The underlying reason is that the dark matter velocity is much larger than the typical sound speed in the fluid, such that the typical energy and momentum transfer for sub-MeV DM cannot match the dispersion relation of a single, on-shell excitation. However, by considering the process of emitting two or more excitations, it is possible to deposit energies larger than meV even with the small momentum transfers characteristic of such light dark matter. The final state excitations are higher-momentum excitations and very nearly back-to-back. The left panel of Fig. 1 illustrates this process.

The purpose of this paper is two-fold. First, we amplify the discussion of Ref. Schutz:2016tid (), providing many more details of the theory utilized for computing the multi-excitation scattering rate. We update the analytic calculation of Ref. Schutz:2016tid () with the measured structure factor for the leading-order scattering rate, and again compare against the available computations of the literature. While neutron scattering data and detailed numerical simulations have been studied in some parts of the multi-excitation phase space, the fluid response for DM with mass below keV rests partially in previously unconsidered regimes of momentum transfer and energy deposition. This therefore requires some theoretical understanding of the rates. Second, we elaborate on the reach for a simplified model of dark matter coupling to nuclei via a new mediator and compare with existing constraints. We additionally consider scattering and absorption via hidden photons, where the final state is a real photon plus a fluid excitation (see right panel of Fig. 1). In this case, the real photon carries away the bulk of the energy, while the fluid excitation absorbs most of the momentum. However, since the net electric charge of a helium atom is screened at the wavelengths of interest, we find that the reach for this case is not competitive with existing stellar constraints.

We introduce the basic elements of the theory for superfluid helium in Sec. II, beginning with a broad introduction to the nature of quasiparticle excitations in the superfluid. In order to calculate the two-excitation process, we employ the correlated basis function formalism, standard in the liquid helium literature, and derive the three-excitation matrix element. App. A provides an alternative formulation in terms of second quantization, and App. B fills in extra details of calculating the three-excitation matrix element. In Sec. III, we turn to a comparison of numerical calculations of the multi-excitation process, applying our results to derive the sensitivity of a liquid helium target to light DM. The results here focus on DM scattering via a mediator that couples to the nucleus. In Sec. IV, we discuss scattering and absorption processes involving a hidden photon, which couples to liquid helium via its polarizability. We conclude in Sec. V.

Ii Theory of superfluid helium

A unique property of helium in the superfluid phase is the nature of the elementary excitations. At long wavelengths (), the elementary degrees of freedom are no longer single-atom excitations. Instead, the elementary excitations are acoustic phonon modes, a collective mode which is equivalent to a density perturbation at long wavelengths. The quasiparticle nature of the phonon modes is also essential for dark matter detection. While phonon modes are present even at temperatures above 2.17 K, the critical temperature for the superfluid phase transition, it is only well below that the width of the phonon mode becomes narrow. In this regime these phonon modes are the only excitations present and they can be thought of as nearly stable quasiparticles. Since a large fraction of the energy deposited in a low-energy dark matter scattering will be in the form of phonons, it is important that these excitations be long-lived states that can propagate to the surface of the liquid and be measured, for a cm volume (or 1 kg) of liquid helium (see Refs. Maris (); WoodsCowley1972 () for a discussion on phonon lifetimes).

In this section, we describe the theory for superfluid helium needed to calculate the production of multiple excitations in the liquid. Due to the strongly-interacting nature of the liquid, the underlying microscopic theory for superfluid helium is not completely understood. However, somewhat phenomenological methods have been proposed which can successfully reproduce many features of the data. The basic idea behind these methods goes back to Feynman in 1954 Feynman54 (), and starts with a posited form for the ground state , or equivalently a wavefunction for an -atom system. While determining the form of the ground state is difficult (though it can be tested by comparison with data), excited states are momentum eigenstates that are written simply as the number density operator acting on the ground state, . This starting point will then allow us to calculate the creation of excitations of the liquid, even without complete knowledge of the full ground state. We will compare this approach with more complete calculations available in the literature. For a thorough discussion of the various theoretical descriptions of excitations in liquid helium, see also Refs. Griffin (); glyde1994excitations (); Feenberg ().

ii.1 Bijl-Feynman relation for single-excitations

Much of our knowledge of the excitations in a strongly-interacting quantum fluid (such as superfluid helium) comes from the dynamic structure function , which describes the response of the liquid to a density perturbation with momentum transfer and energy deposited . For instance, can be directly measured in neutron scattering by measuring the differential cross section:

 d2σdΩdω=b2npfpiS(q,ω), (1)

where and are the initial and final momenta of the scattered neutron, , and is the scattering length of a neutron on an individual helium nucleus.

The dynamic structure function thus depends on the matrix element for the creation of a quasiparticle with momentum and energy . Concretely, is defined as:

 S(q,ω)≡1n0∑β|⟨Ψβ|nq|Ψ0⟩|2δ(ω−ωβ), (2)

where the final states in the scattering are denoted as with energy , the ground state has energy , and . Here is the Fourier transform of the density operator (in real space, ),

 nq≡1√VN∑i=1exp(iq⋅ri), (3)

and are the coordinates of the individual helium atoms in the fluid. We take an arbitrary quantization volume , with the number of He atoms in the volume; physical results will only depend on the average number density . To facilitate some of the later computations, we will occasionally go to the continuum limit by replacing and

The reason neutron scattering (or dark matter scattering) couples to density fluctuations can be understood by considering the potential seen by a neutron in the liquid,

 V(r)=2πbnmn∑iδ(3)(r−ri)=2πbnmnn(r), (4)

assuming a hard-sphere interaction and neutron mass . This is the underlying justification for Eq. (1), which we derive in Sec. III.1 for the case of DM scattering. In particular, we similarly obtain such a potential for dark matter by coupling the DM to helium atoms, with , with and the dark matter mass and scattering length, respectively.

The dynamic structure function is thus crucial to understand the response of superfluid helium to dark matter scattering. While it can be obtained from neutron scattering data at moderate momentum transfer (/Å, corresponding to keV in units where ), a certain level of theoretical control is also possible. As we will see, this theoretical control will be crucial for extrapolating the dynamic structure function to lower momentum transfers, which is necessary to compute the scattering rate when the DM is lighter than keV.

The leading order contribution to is given by the probability to create a single on-shell quasiparticle excitation. One of the earliest theories of the single excitation spectrum, due to Bijl Bijl1940 () and Feynman Feynman54 (), applies the variational method to understand the shape of the dispersion curve. Concretely, the trial wavefunction for a single excitation is given by

 |q⟩=1√n0S(q)nq|Ψ0⟩, (5)

with defined in Eq. (3), and where the static structure function is defined by

 S(q)≡1n0⟨Ψ0|n−qnq|Ψ0⟩. (6)

where is a function only of , and its appearance in the definition of the state ensures that . The left panel of Fig. 2 shows the experimentally measured in helium, which is linear in at small momentum and approaches 1 at high momentum. In the limit of a single excitation, which does not split to multi-excitations, is related to the dynamic structure function by

 S(q,ω) ≈ 1n0∣∣⟨q|nq|Ψ0⟩∣∣2δ(ω−ϵ0(q)) = S(q)δ(ω−ϵ0(q)),

where is the energy of , which we will refer to as the Bijl-Feynman energy.

Since the state in Eq. (5) is by construction orthogonal to the ground state, the variational method dictates that its energy provides an upper bound on the true energy eigenvalue . As we will shown in the next section, the single excitation energy is

 ϵ0(q)≡⟨q|H−E0|q⟩=q22mHeS(q)≥ϵ(q), (8)

with the Hamiltonian and the ground-state energy. The factor comes from the normalization of the states in Eq. (5), and in the limit of a free Bose gas . The Bijl-Feynman theory for produces the single resonance curve shown in the right panel of Fig. 2, and approaches the free-particle quadratic dispersion at high . For comparison, we also show the measured dispersion curve for single-resonance excitations in Fig. 2. We see that the Bijl-Feynman energy agrees roughly with the measured energy at long wavelengths, where the excitations can be identified with sound waves (phonons) with energy , where is the sound speed. In this regime, the static structure factor is then linear in the momentum with

 S(q)≈|q|2mHecs. (9)

However, as the curve reaches a maximum and begins to turn over (the maxon and roton regions), the agreement no longer persists, and is not even qualitatively correct as the dispersion curve reaches a plateau.

The original Bijl-Feynman theory contains, however, no multiphonon response. More generally, the dynamic structure function will contain both the single pole, with strength , and a continuum component, :

 S(q,ω)=Z(q)δ(ω−ϵ(q))+Sm(q,ω), (10)

and the static structure function now satisfies the more general relation . Any large deviation of from indicates that the state defined in Eq. (5) is no longer a good approximation to the single-excitation state, which will be the case in the roton region. The continuum component results from multi-excitation production in the medium. These multi-excitation modes are also important for computing the correct single resonance dispersion curve through radiative corrections to the propagator. It is the multi-excitation response that we will focus on in the rest of this section.

ii.2 Hamiltonian formulation

We now lay out the ingredients to describe phonon interactions, focusing on the elements needed to compute . We follow the correlated basis function formalism, which we briefly review here. This formalism adopts the Bijl-Feynman approach, positing that particle correlations are primarily contained in the ground state wavefunction. Given the exact ground state, excited states are obtained simply with repeated applications of the density operator. Following Ref. JacksonFeenberg (), we define a lowest-order set of basis states using the Bijl-Feynman states:

 |q⟩0 ≡1√n0S(q)nq|Ψ0⟩ (11) |q1,q2⟩0 ≡1√n0S(q1)1√n0S(q2)nq1nq2|Ψ0⟩. (12)

is full ground state of the interacting system. Importantly, the states here are not orthogonal and hence phonon number is not conserved. Instead, the propagating excitations are superpositions of these states. We deal with this complication in the following section.

To compute the energies and matrix elements, we require an interaction Hamiltonian. This Hamiltonian may either be written as the effective theory of a quantum fluid, or in terms of the microscopic degrees of freedom. Let us first consider the fluid Hamiltonian, which directly allows for a second-quantized approach, see e.g. stephen1969raman (). (This discussion most closely follows Ref. Schutz:2016tid ().) Here we elevate the status of the density fluctuation to independent operators which create excitations in the fluid, and consider an effective Hamiltonian for these fluid degrees of freedom,

 H = ∫d3r(12mHev⋅nv+V(n)). (13)

By expanding in the density and velocity fluctuations, the system can be approximated to leading order as a harmonic oscillator with Hamiltonian

 H0 = 12∑qmHen0vq⋅v−q+ϕ(q)nqn−q, (14)

where can be thought of as a momentum dependent force constant. As we show in App. A, this Hamiltonian lends itself to canonical quantization of the variables, and Eq. (14) be can expressed in terms of creation and annihilation operators

 H0=∑qϵ0(q)(a†qaq+12). (15)

The single-excitation energy is simply . The three-excitation interaction vertex and corrected energy eigenvalues can then be obtained by expanding Eq. (13) to higher order in the density and velocity fluctuations. While this setup may be more familiar to a particle physicist, it is less convenient for our purposes. In particular, the ground state in the fluid is nontrivial: in a medium, quantum fluctuations require us to consider an active vacuum, where the asymptotic states of the strongly-interacting fluid are not well-approximated by the free states of a weakly-interacting system. This effect can be accounted for in the second-quantized quantum fluid formalism by correcting the ground state order by order, as we show in App. A, although the calculation is somewhat cumbersome.

In practice, matrix elements are often derived more simply in a first-quantized formulation of the microscopic theory, which has the advantage, as we will see, that knowledge of the ground state is not required to compute the matrix element that we are interested in. Given the energies and vertices computed in this approach, one can of course construct an equivalent second-quantized, quantum fluid Hamiltonian, which may be more convenient for certain scattering and self-energy calculations. The first-quantized microscopic Hamiltonian is given by

 H=∑i(−∇2i2mHe)+V({ri}), (16)

where the sum runs over all particles in the fluid. Writing as the wavefunction corresponding to the ground state , we require that such that is the exact ground state of the full Hamiltonian. For a translationally-invariant system, we thus find that the ground state energy is .

We can show that this formulation also gives the Bijl-Feynman energy in Eq. (8). Using with integration over the coordinates of all atoms, and acting with on the wavefunction ,

 0⟨q|H−E0|q⟩0 =1NS(q)∑i,ℓ∫d3r1...rN ψ0e−iq⋅ri[H−E0](eiq⋅rℓψ0) =12mHe1NS(q)∑i,j,ℓ∫d3r1...rN ψ20(∇je−iq⋅ri)(∇jeiq⋅rℓ)=q22mHeS(q) (17)

where we used and rearranged the derivatives with partial integration. We have also assumed is a properly normalized, real wavefunction, . (Notice that the dependence on the unknown forms of and dropped out.) While we have only computed the average energy for a given state, it can furthermore be shown that approaches an exact eigenstate of in the limit DavisonFeenberg ().

ii.3 Three-excitation vertex

In the previous section, we defined a single excitation state which we regard as an approximately free quasiparticle, as well as multi-excitation states which are products of the single excitations. However, an important subtlety in treating a non-dilute, strongly interacting fluid is that the asymptotic states do not have a well defined particle number. In particular, the states defined so far are not orthogonal, and we must first define an orthogonal basis of states before considering the three-excitation vertex. In other words, to correctly calculate the cross section, we need to compute the matrix element for states which are long-lived compared to the time-scale set by the interaction Hamiltonian. This way the factorization principle allows us to compute the total rate without detailed knowledge of the ultimate fate of the external states in the matrix element.111This is a familiar concept in hard parton scattering in QCD, where we can compute the leading order, total inclusive cross section without detailed knowledge about the shower and hadronization.

A set of orthogonal states can be obtained by performing a Gram-Schmidt rotation to orthogonalize the basis. Concretely, we start with the same single excitation state and then define a set of orthogonal states relative to by

 |q⟩ ≡|q⟩0 (18) |q1,q2⟩ ≡|q1,q2⟩0−∑q′⟨q′|q1,q2⟩0|q′⟩. (19)

In what follows we alway drop the superscript for the single particle state, since it is by construction identical to the corresponding state in the orthogonalized basis. We identify this new basis of states with the orthogonal eigenstates of the quadratic Hamiltonian for the quasiparticles, which will be corrected by the cubic interactions derived below.

The unknown particle correlations of the strongly coupled fluid are now conveniently packaged in the matrix element, which we will discuss later in this section. In the microscopic Hamiltonian of Eq. (16), this overlap term encodes the unknown potential term which dictates the correlations of particles in the ground state. For the quantum fluid effective Hamiltonian in Eq. (13), the same information is encoded in the interactions coming from both the kinetic term, the unknown potential and possible matching terms encoding the unknown short distance physics. (In this sense one may roughly think of the overlap term as a counterterm which enforces the orthogonality of the renormalized states.)

To compute the three-excitation matrix element, we again use with the Hamiltonian given in Eq. (16) and with the ground state energy:

 ⟨q−k,k|δH|q⟩= 0⟨q−k,k|H−E0|q⟩−ϵ0(q) 0⟨q−k,k|q⟩, (20)

In the second term, we have used the leading order energy of the single-excitation state; this three-excitation vertex will itself correct the single-excitation energy at higher order in perturbation theory. The first term in Eq. (20) can be computed directly with the basis states in the previous section:

 0⟨q−k,k|δH|q⟩=1√n30S(q−k)S(k)S(q)∫d3r1...d3rNn∗q−kn∗kψ0(H−E0)nqψ0. (21)

Again, acts on , and after integration by parts plus the fact that satisfies , we can show that

 0⟨q−k,k|δH|q⟩=∑j1√n30S(q−k)S(k)S(q)∫d3r1...d3rN(ψ0)22mHe∇j(n∗q−kn∗k)∇j(nq) (22) =∑j1N√n0S(q−k)S(k)S(q)∫d3r1...d3rN(ψ0)22mHe(−i(q−k)e−i(q−k)⋅rjn∗k−i(k)e−ik⋅rjn∗q−k)(iqeiq⋅rj).

We rewrite the terms above in terms of the static structure function,

 1√N∑i⟨Ψ0|e−iqrinq|Ψ0⟩=1√n0⟨Ψ0|n∗qnq|Ψ0⟩=√n0S(q). (23)

Using this result, we obtain

 0⟨q−k,k|H−E0|q⟩ = q⋅(q−k)S(k)+q⋅kS(q−k)2mHe√N√S(q−k)S(k)S(q) (24)

Next, to directly compute the overlap matrix element requires some working assumption for the form of the ground state wavefunction. Alternatively, one may estimate for this overlap term with a more indirect method. The simplest ansatz which yields the correct long-wavelength behavior and satisfies a certain set of consistency conditions is known as the “convolution approximation.” With this ansatz, one finds JacksonFeenberg (); Feenberg ()

 0⟨q−k,k|q⟩=√S(q−k)S(k)S(q)√N, (25)

which we derive in detail in App. B. It has been shown that using this form gives good agreement with experimental data on neutron scattering. Various improvements to the convolution approximation have been considered (see e.g. PhysRevB.13.3779 ()), though for our approximate, analytic treatment we choose to keep the simplest possibility. This has the main advantage that the formulae of the final answer are very manageable. In particular, utilizing Eq. (20), the full matrix element is then given by

 ⟨q−k,k|δH|q⟩ =q⋅(q−k)S(k)+q⋅kS(q−k)−q2S(k)S(q−k)2mHe√N√S(q−k)S(k)S(q) (26)

Having obtained the three-excitation matrix element, it is now possible to systematically compute the single excitation energy as a perturbation series in this matrix element. To leading order in Brillouin-Wigner perturbation theory brillouinwignerbook (), the eigenstates of the Hamiltonian are

 |Ψq⟩ = |q⟩+12∑p,k|p,k⟩⟨p,k|δH|q⟩ϵ(q)−ϵ0(k)−ϵ0(p)δp+k,q (27) |Ψk,q⟩ = |k,q⟩+12∑p|p⟩⟨p|δH|k,q⟩ϵ(k)+ϵ(q)−ϵ0(p)δp,k+q. (28)

Similarly, the energy of is then given by the recursive relation

 ϵ(q)=⟨Ψq|H|Ψq⟩=ϵ0(q)+12∫d3k(2π)3V|⟨q−k,k|δH|q⟩|2ϵ(q)−ϵ0(q−k)−ϵ0(k). (29)

where we took the continuum limit. To solve for the resummed energy to lowest order, is replaced by inside the integral above. By inserting Eq. (28) in Eq. (2), one can compute the two-excitation contribution to the dynamic structure function to leading order

 Sm(q,ω)=S(q)2∫d3k(2π)3V|⟨q−k,k|δH|q⟩|2(ϵ0(q)−ω)2 δ(ω−ϵ0(k)−ϵ0(q−k)). (30)

This is shown diagrammatically in Fig. 3. Whenever we use this approximation, we use the Bijl-Feynman dispersion relation and the measured form of , both shown in Fig. 2. While this form is enough to obtain a rough estimate of the scattering rate, it clearly has the incorrect structure as it only uses the lowest-order energies.

This deficiency is addressed in many detailed calculations of found in the literature manousakis1986theoretical (); andersenStirling94 (); Krotscheck2015 (), using different approximations for the three-excitation vertex and in defining the multi-excitation states. We note that the approach presented here is not entirely unique in giving reasonable agreement with the data. Recently, a fully self-consistent calculation which resums the corrections due to the three-phonon vertex and gives good agreement with experimental data has been presented by Campbell, Krotscheck and Lichtenegger in Ref. Krotscheck2015 () (hereafter, CKL15). Rather than model the effect of the interactions with a heuristic ansatz for the overlap term, they explicitly include the leading term from the potential in Eq. (16). Operationally, they obtain by recursively solving for the self-energy , which satisfies

 Σ(q,ω)=ϵ0(q)+12∫d3k(2π)3V|⟨q−k,k|δH|q⟩|2ω−Σ(q−k,ω−ϵ0(k))−Σ(k,ω−ϵ0(q−k)). (31)

Using this self-energy, the renormalized energies then match the observed single-excitation energies, and the dynamic structure factor is given by the optical theorem

 S(q,ω)=−1πS(k)ImΣ(q,ω)(ω−ϵ0(q))2+(ImΣ(q,ω))2. (32)

The result for is shown in Fig. 4, which includes both the single and multi-excitation response. We emphasize that the method of Ref. Krotscheck2015 () includes multi-excitation production beyond just the leading order two-excitation production, with the limitation that the multi-excitation production still relies on the three-excitation vertex (in general, higher-point vertices are present). A detailed comparison of this theoretical calculation with inelastic neutron scattering data can be found in Ref. Beauvois2016 (). Accounting for neutrons that scatter multiple times in the liquid, the data is in reasonably good agreement with theory for the multi-excitation component.

As we will discuss in the following section, the results shown in Fig. 4 are in broad agreement with the lowest order calculation of using Eq. (30), although there are significant differences in detailed structure. Where available, we will therefore use the numerical results of CKL15 to compute DM scattering, and use the lowest order results only as a guide to extending CKL15 to low momentum transfer.

Iii Reach for dark matter scattering

We now turn to DM detection with an idealized liquid helium detector, applying our knowledge of the dynamic structure function derived in the previous section. A possible concept for this detector has been shown in mckinseyslac (): the basic idea is that a scattering event creates quasiparticle excitations, which can propagate to the surface of the liquid. At the liquid-gas interface, the quasiparticle has a high probability to eject a free helium atom via quantum evaporation, followed by calorimetric detection of the helium atom. Furthermore, the evaporation process may give a natural amplification technique (with amplification factors of 10), and in principle could be applied for single quasiparticle energies as low as meV.

In this section we use the various results for the dynamic structure function to obtain the rate for DM scattering. We discuss the derivation of the rate given in Schutz:2016tid () in greater detail, considering the expanded calculation of . As a benchmark, we will consider a background-free kg-year exposure. For multi-excitation final states, we take a minimum energy of meV and energies up to 8.6 meV. This upper value on coincides with the upper cutoff of the numerical results we take from CKL15; furthermore, this energy range constitutes the bulk of the response, and the rate falls off rapidly at higher .

The results of this section are applicable to models of dark matter interacting coherently with helium atoms via a new mediator, where we consider both the heavy mediator and light mediator limits. In contrast, in the long wavelength limit the helium atom does not have a net charge for a mediator such as a hidden photon. We will discuss signals related to the hidden photon in Sec. IV.

iii.1 Preliminaries

The total DM scattering rate per unit target mass is given by

 dRdω =1ρHeρXmX∫d3vf(v)∫pi+pf|pi−pf|dq dΓdωdq, (33)

where is the differential scattering rate per incoming DM particle. We denote the initial momentum of the DM and the final momentum , with

 pf=mX√v2−2ωmX  ,  pi=mX|v|. (34)

For the velocity distribution for the dark matter, we assume the standard halo model Maxwellian distribution, boosted to the earth’s frame:

 f(v) =1N(v0,vesc)exp[−(v+ve)2v20]  Θ(vesc−|v+ve|), (35) N(v0,vesc)=π3/2v30[erf(vescv0)−2vescv0exp(−(vescv0)2)] (36)

where km/s, the escape velocity km/s, and we take the average earth’s velocity to be km/s. The normalization factor accounts for the hard cutoff in the distribution at . For the local dark matter density, we take GeV/cm. (Note that this velocity distribution differs somewhat from that used in Ref. Schutz:2016tid (), with more weight at higher initial velocities. This leads to a factor of few larger scattering rate.)

Analogous to the case for neutron scattering Eq. (4), DM in superfluid helium sees the potential

 V(r)=2πbXmX∑iδ(3)(r−ri)=2πbXmXn(r), (37)

where is the DM-helium scattering length. This prescription works both for a light mediator and for a contact operator, where in the former case is momentum dependent. We can then compute the scattering rate with Fermi’s golden rule,

 Γ=2π(2πbXmX)2∫d3pf(2π)3∑β|⟨Ψβ|nq|Ψ0⟩|2δ(Ei−Ef−ωβ). (38)

With a suitable change of variables, the differential rate is then

 dΓdqdω=12n0σX(q)mXqpiS(q,ω), (39)

where we used the definition of in Eq. (2). (The derivation of the neutron scattering rate in Eq. (1) is completely analogous.) The DM-nucleus scattering cross section , where is the DM scattering length. Assuming the DM-nucleus interaction is mediated by a new force carrier , we can express this as

 σX(q)≡⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩σp(fpZ+fn(A−Z))2f2p,mϕ≫q(massive % mediator)σpq4refq4(fpZ+fn(A−Z))2f2p,mϕ≪q(massless mediator), (40)

where we consider the massive and massless mediator limits, and is the DM-proton cross section at a reference momentum transfer . In what follows we take . The expression for DM scattering rate is then

 dRdω=ρX2mHem2X∫d3vf(v)∫pi+pf|pi−pf|dqqpiσX(q)S(q,ω). (41)

iii.2 Scattering rate and reach

Since a full, self-consistent calculation of has been made available in CKL15, we would like to use these results. However, for scattering of light dark matter, the kinematic regime is somewhat different from that of neutron scattering measurements ( keV) and existing simulation data ( eV). In particular, for dark matter in the keV to MeV range, we expect typical momentum transfer and energy deposits given by

 eV≲|q|≲keVandmeV≲ω≲eV. (42)

This is partially outside the regime that was considered in CKL15 and is shown in Fig. 4, which includes 100 eV - 4 keV and meV. The reason for the kinematic mismatch between dark matter and the data is the relatively large velocity of the dark matter compared to the speed of sound in helium, which pushes the interaction away from the linear dispersion phonon regime.

For the time being, we must therefore rely on a theoretically sensible extrapolation to compute the rate for lighter DM. The numerical data in particular shows a scaling in the low region, which we can exploit to extrapolate to lower momenta. This power law can be understood analytically using our approximate expression for in Eq. (30). In the long-wavelength limit ( keV) and at deposited energies   meV, we can take the where and are the momenta of the final state phonons. The matrix element in Eq. (26) then simplifies to

 ⟨q−k,k|δH|q⟩ ≈12mHe√Nq2√S(q)(1−S(k)). (43)

Inserting this in Eq. (30) gives

 S(q,ω)≈116π2q4n0m2Heω2∑i~k2i(1−S(~ki))2, (44)

where the are the solutions to . We show the -dependence of the numerical data from CKL15 in Fig. 5, along with the extrapolation to lower with the power law. For comparison, we also show our own numerical calculations of the leading order using Eq. (30), where we took the Bijl-Feynman dispersion relation and measured form of , each shown in Fig. 2. (While the Bijl-Feynman dispersion relation is strictly speaking not correct for high momenta, we use it to roughly estimate the contribution from the response above 2 meV, as seen in Fig. 4.) In both cases, we see the low- behavior is very well described by a power law.

To indicate the relative importance of this extrapolation for dark matter scattering, we show the values that are most relevant for the DM scattering rate in Fig. 6, compared to the momenta covered by the CKL15 results. What is shown in the average , weighted by the relevant factors in Eq. (41), or more explicitly

 ⟨q⟩≡∫pi+pf|pi−pf|dqq2σX(q)S(q,ω)/∫pi+pf|pi−pf|dqqσX(q)S(q,ω) (45)

using the extrapolated below eV with the power law. Thus the DM rate computed here relies heavily on the extrapolation for DM masses below 50-100 keV, and a dedicated simulation along the lines of CKL15 will eventually be needed in this part of parameter space.

Fig. 7 and Fig. 8 show the spectrum for scattering via a massive and massless mediator, respectively. In both cases, we compare the result using from CKL15 and that using Eq. (30). When computing from Eq. (30), we use the Bijl-Feynman dispersion for excitations, along with the measured ; since this method gives roton/maxon energies which are too high compared to the measured , the structure here is shifted to higher . These differences illustrate the importance of obtaining the correct energies (and widths) of the rotons and maxons, since the rate is clearly dominated by pair-production of these excitations.

The projected best-case sensitivity for DM scattering is shown in Fig. 9, for 1 kg-year exposure and assuming zero background events. The results for both computations of are similar once the rate is integrated over the energy range meV, despite the significant differences in the spectrum. In the same plots, we show the reach if only regular nuclear recoils can be observed down to meV (gray line). (Below meV, we know that the only modes available are quasiparticle (phonon or roton/maxon) modes – see Fig. 2.) In our estimates, we did not include possible backgrounds from scattering of solar neutrinos (see for example Ref. Hochberg:2015fth ()) and coherent photon scattering 2016arXiv161007656R (), which are small for these exposures.

Note our results are consistent with the reach computed in Ref. Schutz:2016tid (), where the results utilizing the CKL15 match exactly (up to the different velocity distributions used). Ref. Schutz:2016tid () also calculated the multi-excitation rate in the leading order approximation, but using a different form of This assumption made it tractable to obtain an analytic result for the rate, but does not include the peaked spectrum from the rotons that we see in Fig. 7 and Fig. 8. However, accounting for a missing symmetry factor of in the analytic results of Ref. Schutz:2016tid () and the different velocity distributions, the reach is similar.

In Fig. 9, we also show contours in for various model-dependent coupling and mass parameters. In particular, the cross section for DM scattering off a single proton or nucleon can be written in the massive and massless mediator limits as

 σp=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩4αXg2nμ2nXm4ϕ,mϕ≫qref(massive mediator)4αXg2nμ2nXq4ref,mϕ≪qref(massless mediator), (46)

for fixed momentum transfer . Here we have written the mediator coupling to the DM and nucleons as and , respectively. (To relate results with the form of the scattering potential given in Eq. (37), we take with , in the limit of .)

Setting aside the cosmological production mechanism for the DM, there are a number of model-dependent existing constraints on light dark matter, in particular for the case of a light mediator. The DM-mediator coupling is bounded from DM self-interactions, which can affect DM halo shapes and small-scale structure. The momentum-transfer weighted self-interaction cross section is given by Tulin:2012wi (),

 σT≈⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩4πα2Xm2Xm4ϕ,mϕ≫mXv(massive mediator)16πα2Xm2Xv4lnmXv22mϕαX,mϕ≪mXv(massless mediator), (47)

where is the velocity of the DM and in the above we have assumed , always valid here. A comparison of observed structure with simulations that incorporate DM self-interactions leads to upper bounds in the ball park of , depending on the system. In particular, observations of dwarf galaxies (with ) allow cross sections as high as Fry:2015rta (); Kaplinghat:2015aga () and comparison with shapes of elliptical galaxies (with gives an upper bound of about Peter:2012jh (). However, we emphasize these bounds can vary by up to an order of magnitude depending on the detailed modeling of structure formation. Furthermore, existing simulations have focused on hard-sphere scattering, and the bounds may be modified significantly in the massless mediator case, where the scattering is dominantly in the forward direction, leading to less isotropization than the hard-sphere scattering case for a given interaction cross-section. Nevertheless, we can use these results to obtain an approximate bound on the DM coupling. Taking and setting , we find

 (48)

where in the light mediator limit we took . Note also that we have assumed , such that quantum mechanical resonance effects can be neglected Tulin:2012wi ().

A new mediator which couples to nucleons is also strongly constrained, for instance by measurements of neutron-nucleus scattering or from stellar cooling constraints. For massive mediators at the MeV scale or heavier, neutron-lead scattering experiments set a constraint of Barbieri:1975xy ()

 gn≲2×10−5(mϕMeV)2. (49)

Combining this with the self-interaction constraint in Eq. (48) gives for an upper limit of , which is well above the cross sections considered here. For reference, we show the cross section for several parameter choices satisfying the self-interaction and neutron-lead scattering constraints in the top panel of Fig. 9.

For the light mediator case, there are also strong constraints from energy loss in helium burning stars, which require for mediator mass  Raffelt:1996wa (). Combined with the self-interaction constraint given in Eq. (48), this would put a strong upper bound on the allowed . However, in both cases the limits are model-dependent and may be uncertain. With this caveat, in Fig. 9, we show for couplings that are roughly consistent with stellar cooling and self-interactions, where for we include the strong -dependence of the bound in Eq. (48).

Additionally, we expect that for the cross sections shown here, the effect of DM stopping in the earth can be neglected (see e.g. Ref. Kavanagh:2016pyr () for a recent detailed analysis of this effect). Assuming an average density of with a chemical composition of Fe (), Si (), O (), Mg () and S (), we estimate the mean free path for DM scattering in the earth to be roughly 9000 km for . The mean free path is therefore larger than the radius of the earth for all cross sections we consider. Moreover, given that every scattering event would only result in a relatively small energy loss, dark matter stopping in the earth can be safely neglected for the cross sections of interest.

Iv Hidden photon processes

A hidden photon is a well-motivated ingredient of many dark matter models, either as a component of the dark matter itself (e.g. Pospelov:2008jk (); Nelson:2011sf (); Arias:2012az ()) or as a mediator for DM interactions. (For a recent review, see Ref. Alexander:2016aln ().) The hidden photon couples to standard model fields through the kinetic mixing operator,

 L⊃κ2FμνF′μν (50)

where is the kinetic mixing parameter and () is the photon (hidden photon) field strength. For a massive hidden photon, this mixing leads to a coupling of the hidden photon with the regular electromagnetic current, , after performing a field redefinition . Here, we consider two scenarios: in the first case, a fermionic DM candidate with keV-MeV mass scatters via the hidden photon mediator. Since our expressions only depend on , with the DM coupling to the hidden photon, for simplicity we set and quote our results only in terms of . In addition, we will calculate absorption of sub-eV mass hidden photons in helium, assuming that the hidden photons constitute the dark matter.

Since the electric charge of the helium atom is screened at long wavelengths (or small momentum transfers keV), the analysis of the previous section no longer applies. Instead, a hidden photon (or photon) couples to the medium by inducing a dipole moment, where the strength of the dipole is determined by the atomic polarizability . Note also that this implies there is negligible difference between the in-medium kinetic mixing and the vacuum kinetic mixing, in contrast to other low-threshold targets like superconductors where in-medium effects substantially affect the rate Hochberg:2015fth ().

Our treatment of this coupling via the polarizability will closely follow Ref. Fetter1972 (), which considered photon scattering in liquid helium. First, we obtain the photon coupling with the medium. To leading order, the target medium is treated as a linear dielectric, with an atomic polarizability (see e.g. Fetter1972 ()) for helium. The polarization of the medium is given by

 P(r)=αn(r)E(r), (51)

where is the number density at helium atoms and is the total electric field in the medium. The interaction Hamiltonian of the polarization with a radiation field is then

 HI=−12∫d3rP(r)⋅Eγ(r). (52)

If the polarization is solely induced by the incident radiation field, then . From the coupling to the number density , this interaction allows for photon scattering by creation of excitations in the liquid. When just a single excitation is emitted, this process is known as Brioullin scattering.222Another possibility is Raman scattering, where in addition to the final state photon, two back-to-back, high momentum phonons are being emitted. However, the rate for Raman scattering is proportional to and is generally three to four orders of magnitude weaker than Brillouin scattering. We neglect it here, but refer to Ref. stephenreview () for a review of both Brillouin and Raman scattering in superfluid helium. Since the sound speed is much smaller than the speed of light, here the phonon excitation only carries a small fraction of the energy, such that the frequency shift in the outgoing photon is minimal.

To obtain the coupling for the hidden photon field, we perform the field redefinition , which gives

 HI=−κα∫d3rn(r)E(r)⋅E′(r) (53)

and is the hidden photon field. From this, we see that the hidden photon couples to a photon and the density field. A DM scattering (or absorption) would thus give rise to both an observable photon and quasiparticle excitation, as shown in the left panel of Fig. 10. The physical interpretation is as follows: an incoming hidden photon must first induce a polarization in the medium, which subsequently relaxes back to the ground state by emitting a photon and a phonon. We calculate the rate for these processes in Sec. IV.1.

Additionally, the polarization vector may be present already if the experimental setup includes a strong external electric field applied in the liquid. In particular, in neutron EDM experiments, superfluid helium is used for storage of the cold neutrons, and a strong electric field is applied to study the neutron spin precession. Recently, a stable electric field as high as kV/cm has been demonstrated Ito:2015hwa (). The interaction Hamiltonian in this case is

 HI=−κα2∫d3rn(r)E0⋅E′(r), (54)

where the external field allows for conversion of the hidden photon into a density perturbation. In this case, there is no final state photon produced, but the kinematics of light DM scattering requires us to consider the multi-excitation final state, analogous to the discussion in previous sections. This process is shown in the right panel of Fig. 10, and we calculate the corresponding rates in Sec. IV.2.

iv.1 Scattering and absorption without an external E-field

We first consider DM scattering in the absence of any external fields. The process is shown in Fig. 10, which also defines our conventions for the kinematic variables. For non-relativistic DM, a typical scattering is characterized by a small deposited energy, but a relatively sizable momentum transfer (). Since the speed of sound in the superfluid is much smaller than the speed of light, nearly all the deposited energy will be carried away by the photon, while the phonon will absorb the momentum:

 ω1≈ω,k2≈qandω2≈0. (55)

To calculate the matrix element, we quantize the electric field of the photon in an arbitrary volume by

 E(r)=i√2V∑<