New constraints and discovery potential of sub-GeV dark matter with xenon detectors

TL;DR

This paper demonstrates that xenon detectors like LUX and LZ significantly improve constraints on sub-GeV dark matter detection, offering enhanced sensitivity and the potential for mass reconstruction of light DM particles.

Contribution

The study introduces new constraints from LUX and LZ on sub-GeV dark matter, surpassing previous limits and highlighting their discovery potential.

Findings

LUX constraints are three orders of magnitude more stringent than previous xenon experiments.

LUX and LZ can distinguish signal from background using photon emission, enabling unambiguous detection.

LZ can reconstruct DM mass with 20% accuracy for particles below 0.5 GeV.

Abstract

Existing xenon dark matter (DM) direct detection experiments can probe the DM-nucleon interaction of DM with a sub-GeV mass through a search for photon emission from the recoiling xenon atom. We show that LUX's constraints on sub-GeV DM, which utilise the scintillation (S1) and ionisation (S2) signals, are approximately three orders of magnitude more stringent than previous xenon constraints in this mass range, derived from the XENON10 and XENON100 S2-only searches. The new LUX constraints provide the most stringent direct detection constraints for DM particles with a mass below 0.5 GeV. In addition, the photon emission signal in LUX and its successor LZ maintain the discrimination between background and signal events so that an unambiguous discovery of sub-GeV DM is possible. We show that LZ has the potential to reconstruct the DM mass with 20% accuracy for particles lighter than 0.5 GeV.

Equations

$$\frac{dR}{d\omega }=\frac{{\rho }_{\mathrm{DM}}}{m_{\mathrm{T}}{m}_{\mathrm{DM}}}{\int }_{v_{\min }}{d}^{3}vvf(\mathbf{v}+{\mathbf{v}}_{\mathrm{E}}){\int }_{E_{\mathrm{R}}^{-}}^{E_{\mathrm{R}}^{+}}d{E}_{\mathrm{R}}\frac{d^2\sigma }{d\omega d{E}_{\mathrm{R}}},$$

$$\frac{d^2\sigma }{d\omega d{E}_{\mathrm{R}}}=\frac{4\alpha }{3\pi }\frac{{\left|f(\omega )\right|}^2}{\omega }\frac{E_{\mathrm{R}}}{m_{\mathrm{T}}}\frac{d\sigma }{d{E}_{\mathrm{R}}},$$

$$\frac{d\sigma }{d{E}_{\mathrm{R}}}=\frac{m_{\mathrm{T}}{\sigma }_{\mathrm{SI}}^0}{2{\mu }_n^2{v}^2}{F}_{\mathrm{T}}^{\mathrm{SI}}(E_{\mathrm{R}})F_{\mathrm{med}}(E_{\mathrm{R}}),$$

$$\frac{d^{2}R}{dS1dS2}=ϵ(S1,S2)\int d\omega \frac{dR}{d\omega }\mathrm{pdf}(S1,S2|\omega ),$$

$$n_q$$

$$=n_{\mathrm{ion}}+n_{\mathrm{ex}},$$

$$n_e=n_{\mathrm{ion}}(1-r),$$

$$n_{\gamma }$$

$$=n_{\mathrm{ion}}(r+\alpha )$$

$$r=\frac{L_y-\alpha Q_y}{L_y+Q_y}.$$

$$r=1-\frac{\log (1+\xi n_{\mathrm{ion}})}{\xi n_{\mathrm{ion}}},$$

$$Q_y=\frac{1}{\xi E}\log \left(1+\frac{\xi E}{W(1+\alpha )}\right).$$

Full text

††preprint: KCL-PH-TH/2017-16

I Introduction

Identifying the nature of dark matter (DM) remains one of the most compelling problems in astroparticle physics. Motivated by the weakly interacting massive particle (WIMP) paradigm, DM direct detection experiments have traditionally concentrated on the 5 GeV to 10 TeV mass range. The non-detection of DM in this range has led to significant theoretical efforts focussing on lighter particles. Initially, sub-GeV DM was mainly considered in the context of the 511 keV gamma-ray anomaly observed by INTEGRAL Boehm and Fayet (2004); Boehm et al. (2004a, b), but more recently, it has been recognised that sub-GeV DM is generic in many other scenarios, e.g. Pospelov et al. (2008a); Hall et al. (2010); Hochberg et al. (2014); D’Agnolo and Ruderman (2015); Pappadopulo et al. (2016); D’Agnolo et al. (2017).

This theoretical activity has motivated both new direct detection experiments for sub-GeV DM Essig et al. (2016); Hochberg et al. (2016a, 2017a); Derenzo et al. (2017); Hochberg et al. (2016b, 2017b); Bunting et al. (2017); Guo and McKinsey (2013); Schutz and Zurek (2016); Knapen et al. (2017); Essig et al. (2017), and new searches with existing experiments Pospelov et al. (2008b); Essig et al. (2012a, b); An et al. (2015); Lee et al. (2015); Bloch et al. (2017). The major obstacle faced by low mass searches is that the energy deposited in a detector by sub-GeV DM is small. For instance, the maximum recoil energy imparted by DM (of mass $m_{\mathrm{DM}}$) to a nucleus (with mass number $A$) in elastic scattering is $E^{\rm{max}}_{\rm{R}}\approx 0.1~{}\mathrm{keV}\cdot(131/A)\,(m_{\mathrm{DM}}/1~{}\mathrm{GeV})^{2}$. The nuclear recoil energy threshold of dual-phase xenon detectors $(A_{\rm{Xe}}\simeq 131)$ is approximately $1$ keV, implying that they are limited to $m_{\mathrm{DM}}\gtrsim 3$ GeV. The lighter nuclei and lower energy thresholds employed in the CRESST Angloher et al. (2016, 2015), DAMIC Aguilar-Arevalo et al. (2016), EDELWEISS Hehn et al. (2016) and SuperCDMS Agnese et al. (2016, 2017) detectors allow them to probe lower masses, with current exclusion limits reaching down to $m_{\mathrm{DM}}\simeq 0.5$ GeV. Unfortunately, the push to a lower energy threshold often comes with the loss of good discrimination between background and DM events, limiting their ability to make an unambiguous discovery of sub-GeV DM.

Reference Kouvaris and Pradler (2017) demonstrated that existing xenon detectors can probe sub-GeV DM through a new signal channel: a search for the irreducible photon emission from a polarised xenon atom, caused by the displacement of the nucleus and electron charges after the xenon nucleus recoils, and derived constraints from the XENON10 and XENON100 S2-only searches Angle et al. (2011); Aprile et al. (2016a). In this $\mathrm{DM}+\mathrm{Xe}\to\mathrm{DM}+\mathrm{Xe}+\gamma$ inelastic scattering process with a photon in the final state, the maximum photon energy is $\omega^{\rm{max}}\approx 3~{}\mathrm{keV}\cdot(m_{\mathrm{DM}}/1~{}\mathrm{GeV})$. In this paper, we show that the LUX dual-phase xenon detector can also probe sub-GeV DM with the more powerful $\mathrm{S1}+\mathrm{S2}$ search, which is sensitive to photon energies $\omega\gtrsim 0.3$ keV.111The S1 and S2 signals are defined more carefully in Sec. III. We calculate the parameter space excluded with data from LUX’s two WIMP searches (WS2013 Akerib et al. (2016a) and WS2014-16 Akerib et al. (2017a)) and show that the LUX constraints are up to three orders of magnitude more sensitive than the XENON10 and XENON100 S2-only searches considered in Ref. Kouvaris and Pradler (2017). This is because the photon’s energy is large enough to produce detectable scintillation and ionisation charge signals, with the result that events from the fiducial volume, where the background is lower Marrodan Undagoitia and Rauch (2016), can be selected. Moreover, the good discrimination between background and signal events based on the scintillation and ionisation signals is retained. This further reduces the background rate and importantly, allows for an unambiguous detection of DM to be made. We demonstrate this explicitly for LZ Akerib et al. (2015), where we calculate its sensitivity and show that an experiment under construction has the potential to accurately reconstruct the parameters of sub-GeV DM.

II Photon emission scattering rate

The differential rate for a DM particle to undergo two-to-three scattering with a nucleus of mass $m_{\mathrm{T}}$ is

[TABLE]

where $\rho_{\rm{DM}}=0.3~{}\mathrm{GeV}/{\mathrm{cm}}^{3}$ is the local DM density and $f(\mathbf{v})$ is the DM velocity distribution in the Galactic frame, which we assume is a Maxwell-Boltzmann (MB) distribution with a cut-off at $v_{\rm{esc}}=544~{}\mathrm{km}/\mathrm{s}$ and most probable speed of $v_{0}=220~{}\mathrm{km}/\mathrm{s}$. We boost from the Galactic to the Earth reference frame with $\mathbf{v}_{\rm{E}}$ McCabe (2014); Lee et al. (2013). Small deviations from a MB distribution are likely, as seen in numerical simulations, e.g. Bozorgnia et al. (2016); Kelso et al. (2016); Sloane et al. (2016), and predicted by Earth-scattering effects, e.g. Collar and Avignone III (1992, 1993); Kavanagh et al. (2017), but we do not consider them in this work. The limits of integration are found from three-body kinematics, $E_{\rm{R}}^{\pm}=\mu_{\mathrm{T}}^{2}v^{2}/m_{\mathrm{T}}\cdot\left[1-v_{\rm{min}}^{2}/(2v^{2})\pm\sqrt{1-v_{\rm{min}}^{2}/v^{2}}\right]$, where $v_{\mathrm{min}}=\sqrt{2\omega/\mu_{\mathrm{T}}}$ is the minimum DM speed required for a photon to have energy $\omega$, and $\mu_{\mathrm{T}}$ is the DM-nucleus reduced mass.

The photon emission cross-section (derived in Ref. Kouvaris and Pradler (2017)) is

[TABLE]

where $\alpha$ is the fine-structure constant, $f(\omega)=f_{1}(\omega)+if_{2}(\omega)$ are atomic form factors ato , and $d\sigma/dE_{\rm{R}}$ is the DM-nucleus cross-section for elastic scattering. The price to pay for photon emission is a factor $E_{\rm{R}}/m_{\mathrm{T}}$, resulting in a $\mathcal{O}(0.1~{}\mathrm{keV}/100~{}\mathrm{GeV})\simeq\mathcal{O}(10^{-9})$ suppression factor.

We parameterise the DM-nucleus cross-section as

[TABLE]

where $\mu_{n}$ is the DM-nucleon reduced mass and $F^{\rm{SI}}_{\rm{T}}$ is the nuclear form factor. It is an excellent approximation to evaluate $F^{\rm{SI}}_{\rm{T}}$ at $E_{\rm{R}}=0$ keV. We focus on spin-independent (SI) interaction with equal interaction strength with protons and neutrons so $F^{\rm{SI}}_{\rm{T}}=A^{2}$ Vietze et al. (2015). Finally, $F_{\rm{med}}(E_{\rm{R}})$ is a factor that depends on the mass of the particle mediating the interaction. In the heavy mediator limit, $m_{\rm{med}}\gg q$, where $q=\sqrt{2m_{\mathrm{T}}E_{\mathrm{R}}}\sim 3~{}\mathrm{MeV}\cdot(m_{\mathrm{DM}}/1~{}\mathrm{GeV})$ is the momentum transfer, $F_{\rm{med}}=1$ and hence, $\sigma^{0}_{\rm{SI}}$ is the usual DM-nucleon cross-section that is constrained in SI analyses. In the light mediator limit, $m_{\rm{med}}\ll q$, $F_{\rm{med}}=q_{\rm{ref}}^{4}/q^{4}$. In this case, $\sigma^{0}_{\rm{SI}}$ must be defined at a reference value of $q$; we take $q_{\rm{ref}}=1~{}\mathrm{MeV}$, the typical size of $q$ for $m_{\mathrm{DM}}\lesssim 1$ GeV.

III LUX and LZ detector simulations

Dual-phase xenon detectors do not directly measure energy. Rather, they measure the ‘S1’ and ‘S2’ signals, proportional to the initial scintillation and ionisation charge respectively, produced by an energy deposition Chepel and Araujo (2013).

The event rate in terms of the observable signals is

[TABLE]

where $\epsilon(\mathrm{S1},\mathrm{S2})$ represent detection efficiencies and we determine $\mathrm{pdf}(\mathrm{S1},\mathrm{S2}|\omega)$ with a Monte Carlo simulation of the detector. Our simulations are based on the Noble Element Simulation Technique (NEST) Szydagis et al. (2011, 2013); Lenardo et al. (2015) and following Ref. Szydagis et al. (2013), we assume that the electron and photon yields, $Q_{y}$ and $L_{y}$ respectively, from keV-energy depositions from beta particles and gamma rays are the same (collectively, electronic recoils). For an energy $E$, the mean S1 and S2 values are related to the yields through $\mathrm{S1}=g_{1}L_{y}E$ and $\mathrm{S2}=g_{2}Q_{y}E$, where $g_{1}$ and $g_{2}$ are proportionality (or ‘gain’) factors. Following theoretical arguments, we assume that the $\mathcal{O}(0.1)~{}\mathrm{keV}$ nuclear recoil associated with the $\mathrm{DM}+\mathrm{Xe}\to\mathrm{DM}+\mathrm{Xe}+\gamma$ inelastic scattering process does not produce an observable signal Sorensen (2015). There are proposals to test this assumption with new low-energy calibration techniques Verbus et al. (2017).

First, we describe our simulation for electronic recoils (ERs), where our input is $Q_{y}$. Above 1.3 keV, we fit $Q_{y}$ to the central values of LUX’s tritium calibration data Akerib et al. (2016b). Below 1.3 keV, we fit to the central values from LUX’s calibration with ${}^{127}\mathrm{Xe}$ xe (1). There are no calibration data below 0.19 keV so we assume that $Q_{y}$ is zero below this energy. We self consistently determine $L_{y}$ through the relation $n_{q}/E=Q_{y}+L_{y}$, where $n_{q}=E/13.7~{}\mathrm{eV}$ is the total number of quanta from an energy $E$. Our $L_{y}$ agrees perfectly with LUX’s $L_{y}$ calibration data above 1.3 keV Akerib et al. (2016b). As with $Q_{y}$, we assume that $L_{y}$ is zero below 0.19 keV. Our yields are also in good agreement with data from the PIXeY xenon detector Bou ; Boulton et al. (2017). We include recombination fluctuations, generating the recombination probability $r$ and the fluctuations as described in Ref. Akerib et al. (2017b) with parameters $\sigma_{p}=0.07$ and $\alpha=0.2$. Additional parameterisations of $Q_{y}$ are investigated in Appendices A and B.

Second, we summarise our simulation for nuclear recoils (NRs). This is used to check that our simulations correctly reproduce published LUX NR results and also to calculate the ${}^{8}\rm{B}$ solar neutrino signal in LZ. As input, we use a $Q_{y}$ and $L_{y}$ parameterisation that fits LUX’s D-D calibration data Akerib et al. (2016c) and assume the Lindhard model with $k=0.174$ Akerib et al. (2016c). We include Penning quenching as in Ref. Akerib et al. (2016c), model the recombination probability following the Thomas-Imel model with parameters in Ref. Bailey (2016), and again use $\sigma_{p}=0.07$ to model recombination fluctuations.

Unless stated otherwise, S1 and S2 refer to position corrected values, where S1 is normalised to the centre of the detector and S2 to the top of the liquid. We take into account the variation of $g_{1}$ with height within the detector (we ignore radial variations) using results in Ref. Dobi (2014) for LUX and projections for LZ in Ref. Akerib et al. (2015). For LUX (LZ), we use an electron lifetime of 800 (3000) ms and an electron drift speed 1.5 mm/$\mu$s for both.

For LUX WS2013, we use the parameter values from Refs. Akerib et al. (2016a, d): $g_{1}=0.117~{}\mathrm{phd}/\gamma$, $g_{2}=12.1~{}\mathrm{phd}/e^{-}$, an extraction efficiency of $49\%$; the S1 detection efficiency from Ref. Bailey (2016); we allow events that satisfy $\mathrm{S1}_{\mathrm{raw}}>1~{}\mathrm{phd}$, $\mathrm{S2}_{\rm{raw}}>165~{}\mathrm{phd}$; and compare against events measured with $\mathrm{radius}<18$ cm to set an exclusion limit.

LUX WS2014-16 was more complicated owing to the changing conditions throughout the run. We do not attempt to model the changing conditions with time or the spatially varying electric drift field. Instead we take a simplified approach and assume the average values from Ref. Akerib et al. (2017a): $g_{1}=0.1~{}\mathrm{phd}/\gamma$, $g_{2}=18.9~{}\mathrm{phd}/e^{-}$, an extraction efficiency of $73\%$; the S1 efficiency from Ref. Bailey (2016) and the S2 efficiency from Ref. szy ; we allow events that satisfy $\mathrm{S1}_{\mathrm{raw}}>1~{}\mathrm{phd}$, $\mathrm{S2}_{\rm{raw}}>200~{}\mathrm{phd}$. We set an exclusion limit with events measured more than 1 cm from the radial fiducial volume boundary.

For LZ, we use parameter values recommended in Ref. Akerib et al. (2015): $g_{2}=50~{}\mathrm{phd}/e^{-}$, an extraction efficiency of $100\%$ and allow events that satisfy $\mathrm{S2}_{\rm{raw}}>250~{}\mathrm{phd}$. The S1 signal is the main determinant of the energy threshold so we show results taking the upper (lower) values of the range in Ref. Akerib et al. (2015), namely $g_{1}=0.1\,(0.05)~{}\mathrm{phd}/\gamma$, and allow events that satisfy $\mathrm{S1}>2\,(3)$ phd.

Before presenting the limits on sub-GeV DM, we demonstrate that our simulations accurately reproduce published LUX results. First, we derive the efficiency as a function of energy for ERs and NRs, shown in the insets of Fig. 1, and compare against the LUX values (black triangles). Good agreement is found; within 5% above 2 keV for both runs. Below 2 keV, we slightly underestimate the published efficiencies, reaching an underestimation of 50% at 1.1 keV. Results for a direct comparison of the ER efficiency are not available so instead, we compare with the tritium calibration run Akerib et al. (2016b), which had only slightly different parameters from WS2013 and WS2014-16. In the calibration run, the ER efficiency was 50% at 1.24 keV, consistent with 1.13 and 1.25 keV that we find for WS2013 and WS2014-16, respectively.

Lastly, in Fig. 1, we compare our ER and NR bands, indicated by the blue and red lines respectively, against the LUX bands, indicated by squares and triangles at 2 phd intervals (we use updated LUX WS2013 bands Sha ). We find good agreement in both the central position and the width of the bands.

IV LUX constraints

We use two methods to derive constraints on the DM-nucleon cross-section from the WS2013 Akerib et al. (2016a, d) and WS2014-16 Akerib et al. (2017a) WIMP searches.

The first method is a cut-and-count (CC) approach, the simplest and most conservative approach that treats all measured events as signal events. For each mass, we calculate the signal region that contains 90% of the DM events that pass all cuts. The black lines in Fig. 1 give examples of this region for different values of $m_{\mathrm{DM}}$ in the heavy mediator limit. For $m_{\mathrm{DM}}=0.1$ GeV, all of the DM events lie far below the ER band, where background events are expected to lie. For this mass, the mean S1 signal is below the S1 threshold so only the tail of the upward S1 fluctuations is measured. The S1 signal can fluctuate upwards in our simulation from the binomial modelling of the initial number of ions and excitons, the binomial modelling of the detection of photons by the photomultiplier tubes, the Gaussian resolution of the detector or finally, from the modelling of recombination fluctuations.222See Ref. Dobi (2014) for an extended discussion on fluctuations in LUX. As $m_{\mathrm{DM}}$ increases, more of the DM contour overlaps the ER band because fewer upward S1 fluctuations are probed. The signal more closely follows the ER band but there is still a small offset. For clarity, we do not show the contours for a light mediator in Fig. 1. They are similar but extend to slightly smaller S1 values e.g. to $\mathrm{S1}=4.4\,(5.9)$ phd for $m_{\mathrm{DM}}=0.6\,(1.0)$ GeV.

We count all of the observed events within the 90% DM contour and use Poisson statistics to set a 90% C.L. exclusion limit. This is shown as the pink and purple lines in Fig. 2 for the WS2013 and WS2014-16 WIMP searches respectively. The kinks around $m_{\mathrm{DM}}\sim 0.3$ GeV occur as measured events suddenly enter the signal region.

The second more powerful method to derive a constraint on the DM-nucleon cross-section uses a profile likelihood ratio (PLR) test. This takes into account the S1 and S2 information of each event and allows for the results from WS2013 and WS2014-16 to be combined. Unlike the CC method, the background signal must be quantified. We adopt a simple model that assumes the ER background rate is flat in energy, while ignoring subdominant contributions from neutrons and ${}^{8}\mathrm{B}$ neutrinos. Additionally, for the WS2013 search we include a component for the decays of ${}^{127}\mathrm{Xe}$, which contributed in WS2013 but not in WS2014-16. This simple model provides a good fit to LUX data Szydagis (2016); Akerib et al. (2017c) (a comparison is provided in Appendix C). The 90% C.L. combined PLR limit is calculated following Ref. Cowan et al. (2011) with an un-binned extended likelihood function Barlow (1990). We follow the safeguard method in Ref. Priel et al. (2017) to minimise the effect of background mismodelling. The amplitudes of the background components in each run are treated as a nuisance parameter. The 90% C.L. limit is shown as the red line in Fig. 2. It is similar to the CC limit at low mass, where the DM signal region is far from the background region. At higher mass, the PLR limit is significantly stronger. At $m_{\mathrm{DM}}=1~{}\mathrm{GeV}$, the limit corresponds to 5.0 signal events in WS2014-16 (for heavy and light mediators).

In Fig. 2, we also show 90% C.L. exclusion limits from CRESST-II Angloher et al. (2016) and XENON100’s S2-only analysis Aprile et al. (2016a), and a constraint from Earth heating by DM annihilation Mack et al. (2007) (which does not apply if DM annihilation does not occur). For CRESST-II, we use Lise’s public data Angloher et al. (2017) and calculate a limit with the maximum gap method Yellin (2002). XENON100 observed a large number of events (13560) so we use a CC analysis with events in the range 80-1000 PE to set a limit. Our XENON100 limit is slightly weaker than in Ref. Kouvaris and Pradler (2017) because we adopt the $Q_{y}$ used in our LUX analysis, which has a cut-off at 0.19 keV. Before this work, CRESST-II and XENON100 (S2-only) provided the most stringent direct detection constraints on the DM-nucleon cross-section in this mass range. The LUX limits are significantly stronger and already reach the principal reach projected in Ref. Kouvaris and Pradler (2017). This is because the displacement of signal and background regions, not previously considered, further reduces the background rate below the value considered in Ref. Kouvaris and Pradler (2017).

V LZ sensitivity projection

The next generation of dual-phase detectors, namely LZ, XENON1T/XENONnT Aprile et al. (2016b) and PandaX-II/PandaX-xT Cao et al. (2014), will be bigger than LUX while having a smaller background rate. We focus on sensitivity projections for LZ because it has the most detailed design and performance studies Akerib et al. (2015); Mount et al. (2017). Within LZ’s fiducial volume, ${}^{8}\mathrm{B}$ (pp and ${}^{7}\mathrm{Be}$) solar neutrinos produce the dominant NR (ER) background events. Figure 3 (upper panel) shows a simulation of the events seen with the more sensitive LZ scenario that we consider ($g_{1}=0.1~{}\mathrm{phd}/\gamma$, $\mathrm{S1}\geq 2~{}\mathrm{phd})$, together with the ER and NR efficiencies. Assuming a 5.6 tonne fiducial mass, 1000 days of data and the neutrino fluxes from Ref. Bergstrom et al. (2016), we find that 49.8 ${}^{8}\mathrm{B}$ events are expected. In the less sensitive scenario ($g_{1}=0.05~{}\mathrm{phd}/\gamma$, $\mathrm{S1}\geq 3~{}\mathrm{phd})$, the ER and NR efficiencies (not shown) are shifted to higher energies. The ER (NR) efficiency is 1% at 0.94 (2.3) keV respectively, so that only 2.7 ${}^{8}\mathrm{B}$ events are expected.

We quantify LZ’s sensitivity by calculating the median cross-section for LZ to make a discovery at $3\sigma$ (or greater) significance. This is shown for the two LZ scenarios by the black and grey lines in Fig. 2. We use a PLR test and include a 2.5% (1%) uncertainty on the ${}^{8}\mathrm{B}$ ($\mathrm{pp}+^{7}\mathrm{Be}$) flux Bergstrom et al. (2016). For both scenarios, these cross-sections correspond to approximately $5\,(15)$ expected signal events at $m_{\mathrm{DM}}=0.3\,(1.0)$ GeV. For the more sensitive scenario, we also show the $\pm 1\sigma\,(\pm 2\sigma)$ containment region in green (yellow). The more sensitive LZ scenario results in a factor 20–50 improvement compared to LUX, while the less sensitive scenario leads to only a small improvement above $m_{\mathrm{DM}}=0.3$ GeV.

For the more sensitive LZ scenario, we investigate the precision with which LZ can reconstruct the parameters of sub-GeV DM. Figure 3 (lower panel) shows examples of two reconstructions where a high-significance ($>5\sigma$) detection of DM is made. All of the signal and background events used in the reconstruction in the upper benchmark ($\bigstar$ DM benchmark) are displayed in the upper panel of Fig. 3. For both benchmarks, the mass is reconstructed with $\simeq 20\%$ accuracy.

VI Summary

Upcoming xenon detectors will provide opportunities to search for signals beyond the standard DM-nucleus interactions. Previous studies have investigated recoils induced from: solar neutrinos Baudis et al. (2014); Cerdeno et al. (2016), supernova neutrinos Chakraborty et al. (2014); Davis (2016); Lang et al. (2016), nuclear DM Hardy et al. (2015); Butcher et al. (2016), products from DM annihilation Cherry et al. (2015); and inelastic nucleus scattering Baudis et al. (2013); McCabe (2016).

We have investigated photon emission from the recoiling atom, another non-standard signal that allows dual-phase xenon detectors to probe sub-GeV DM. We have demonstrated that the LUX constraints are approximately three orders of magnitude more constraining than the S2-only limits from XENON100, and extend to smaller masses than the CRESST-II limit. In addition, a future experiment such as LZ can accurately reconstruct the parameters of sub-GeV DM since dual-phase detectors maintain the discrimination between background and signal events.

Acknowledgements.

C.M. thanks Jelle Aalbers, Gianfranco Bertone, Andrew Brown, Patrick Decowski, Thomas Edwards, Chamkaur Ghag, Achim Gutlein, Federica Petricca, Ludwig Rauch, Florian Reindl, Peter Sorensen and Matthew Szydagis for discussions and correspondence. C.M. gratefully acknowledges support from the Netherlands Organisation for Scientific Research (NWO); and the Science and Technology Facilities Council (STFC) Grant ST/N004663/1.

Appendix A Signal models for electronic recoils

To aid the reproducibility of our results, we here provide the details of our signal generation model for electronic recoil (ER) events. We begin by introducing the general formalism before presenting assumptions specific to the model used in the main part of the paper. We then introduce two additional models based on the Thomas-Imel box model Thomas and Imel (1987). The material discussed below only addresses the mean signal yields. We leave the details of our model for fluctuations to the main part of the paper.

The number of quanta $n_{q}$ for ER events, in terms of the number of photons $n_{\gamma}$ and electrons $n_{e}$, or in terms of the number of ions $n_{\rm{ion}}$ and excitons $n_{\rm{ex}}$, is

[TABLE]

where $n_{q}=E/W$, $E$ is the energy and we take $W=13.7$ eV.

The number of electrons is related to the electron yield $Q_{y}$ (also referred to as the charge yield) by $n_{e}=Q_{y}E$. The starting point for our detector simulation for ERs takes $Q_{y}$ as input. We fit $Q_{y}$ to LUX’s tritium calibration data Akerib et al. (2016b) above 1.3 keV, while below this, we fit to the central values from LUX’s low-energy calibration with ${}^{127}\mathrm{Xe}$ xe (1). These data points are shown by the blue and red data points in the upper panel of Fig. 4, respectively. In addition, the yellow data points and lilac boxes in Fig. 4 show the low-energy calibration data taken with the PIXeY and neriX xenon detectors. The PIXeY data are from the decays of ${}^{37}\mathrm{Ar}$ Bou ; Boulton et al. (2017), while the neriX data were taken at 190 V/cm and we show the dominant systematic uncertainty Goetzke et al. (2016). We do not include the PIXeY or neriX data in our fits but they are consistent with the ${}^{127}\mathrm{Xe}$ data and tritium data.

The number of electrons is related to the number of ions through

[TABLE]

where $r$ is the fraction of ions that undergo recombination. From the sum rule in Eq. (5), we must have that

[TABLE]

where we have introduced the parameter $\alpha=n_{\rm{ex}}/n_{\rm{ion}}$. Values for $\alpha$ used in the literature typically fall between $0.06$ and $0.2$ Chepel and Araujo (2013).

A.1 Model 1: Interpolation through the data (used in the main part of the paper)

For the model used in the main part of the paper, $Q_{y}$ is determined by tracing a line (on a log axis) through the central points of LUX’s tritium and ${}^{127}\mathrm{Xe}$ data points. This is shown by the black solid line in Fig. 4. There is no physical basis behind this model and we include an unphysical cut-off in $Q_{y}$ at the energy of the lowest data point (0.19 keV). The photon yield $L_{y}$ (also referred to as the light yield) is then straightforwardly determined through the relation $1/W=Q_{y}+L_{y}$ [a rearrangement of Eq. (5)]. The solid black line in the second panel from the top of Fig. 4 shows this parameterisation of $L_{y}$. It passes through the central values of $L_{y}$ from LUX’s tritium calibration data, shown by the blue data points. The lowest ${}^{127}\mathrm{Xe}$ data point in the upper panel of Fig. 4 satisfies $Q_{y}\approx 1/W$, which explains why $L_{y}$ tends to zero at approximately 0.19 keV.

From simple algebraic manipulation of Eqs. (5) to (9), we can express the recombination probability $r$ as

[TABLE]

At higher energies, we follow LUX and fix $\alpha=0.2$ Akerib et al. (2017b). Under these assumptions, $r>0$ for $E>1.1$ keV (see the solid line in the recombination panel in Fig. 4). For energies smaller than this, we fix $r=0$ by requiring that $\alpha=L_{y}/Q_{y}$. With this approximation, $\alpha$ smoothly decreases from 0.2 to zero at approximately 0.19 keV (shown by the solid line in the bottom panel of Fig. 4), as it must to ensure that $n_{\gamma}=0$ at approximately 0.19 keV [cf. Eq. (9)].

A.2 Model 2: Thomas-Imel box model

The previous model was ad-hoc in that it was chosen to pass through the central value of the data points without any relation to a physically motivated model. In particular, the cut-off at approximately 0.19 keV is unphysical as it is more reasonable to expect that quanta are produced all the way to energies $\mathcal{O}(W)$. We therefore now explore the implications of a more physically motivated model of recombination: the Thomas-Imel box model Thomas and Imel (1987). In this model

[TABLE]

where $\xi$ is a free parameter.

This model for $r$ allows us to rewrite the charge yield as

[TABLE]

In this case, by fixing $\xi$ and $\alpha$, we determine $Q_{y}$, $r$ and $L_{y}\,(=1/W-Q_{y})$. We perform a $\chi^{2}$ fit to the tritium and ${}^{127}\mathrm{Xe}$ data below 10 keV to fit $\xi$ and $\alpha$, finding $\alpha=0.06$ and $\xi=0.0065$. The resulting values of $Q_{y}$, $L_{y}$, $r$ and $\alpha$ for this model are shown by the dotted lines in Fig. 4. This model slightly underestimates (overestimates) the $Q_{y}$ $(L_{y})$ low-energy data points.

A.3 Model 3: Thomas-Imel inspired model with energy-dependent coefficients

The low-energy ${}^{127}\mathrm{Xe}$ and PIXeY data are consistent with $n_{e}\approx n_{q}$ (or equivalently $Q_{y}\approx 1/W\approx 73~{}e^{-}/\mathrm{keV}$), while in the Thomas-Imel model $Q_{y}\to 1/(W(1+\alpha))$ as $E\to 0$. Therefore to improve the fit with the ${}^{127}\mathrm{Xe}$ and PIXeY data, we require that $\alpha\to 0$ at low energy.

We therefore modify the Thomas-Imel model to make the parameters in the Thomas-Imel model energy dependent at low energy. We define $\alpha=\alpha_{0}\cdot\mathrm{erf}(\alpha_{1}E)$ and $\xi=\xi_{0}\cdot\mathrm{erf}(\xi_{1}E)$. The justification for this parameterisation is simply to provide a way to smoothly force the parameters to zero as $E\to 0$ in order that the fit with the ${}^{127}\mathrm{Xe}$ and PIXeY data is improved. We perform a $\chi^{2}$ fit to obtain $\{\xi_{0},\xi_{1},\alpha_{0},\alpha_{1}\}$. By construction, in this parameterisation $Q_{y}$ and $L_{y}$ now provide a good fit to the low-energy data points. Comparing the dashed and solid lines in Fig. 4, we see that Model 3 is similar to Model 1 but it avoids the cutoff at 0.19 keV in $Q_{y}$ and $L_{y}$, and the sharp transitions in $\alpha$ and $r$ at 1.1 keV.

Finally, we have also shown the $Q_{y}$ parameterisation adopted in the 2017 LZ technical design report (TDR) Mount et al. (2017) by the dot-dashed green line in the upper panel of Fig. 4. Both Model 1 and Model 3 closely resemble this parameterisation.

Appendix B Exclusions limits and discovery potential with different signal models

We now explore the impact of the different signal models on the LUX exclusion limits and the LZ discovery potential. The results for Model 1 are presented in the main part of the paper so we here focus on the results for Model 2 and Model 3.

We first focus on the exclusion limits from LUX WS2013. In the upper panels of Fig. 5, we have recalculated the ER band (blue lines) and 90% signal regions for two values of the dark matter mass ($m_{\mathrm{DM}}$). The NR model is the same across all plots. With Model 2, the ER band extends to lower $\mathrm{S2}$ values for small $\mathrm{S1}$ values compared to Model 1 and Model 3. Meanwhile, the dark matter signal regions show only minor changes. Therefore the discrimination between signal and background is not as high in Model 2. In Models 2 and 3, the ER efficiency extends to lower energies compared to Model 1 (a direct comparison is made in the inset of the upper right panel of Fig. 5) because in these models $L_{y}$ is non-zero below 0.19 keV. The result of the higher efficiency at lower energies is that the LUX limits are more stringent for Models 2 and 3 at lower values of $m_{\mathrm{DM}}$ (see the limit plot in upper right panel of Fig. 5).

We next explore the discovery potential for the more sensitive LZ scenario. In the lower panels of Fig. 5, we have again recalculated the ER band (blue lines) and 90% signal regions for two values of the dark matter mass. We again find that the ER band in Model 2 extends to lower $\mathrm{S2}$ values for small $\mathrm{S1}$ values compared to the other models. The signal regions are again somewhat similar, although the effect of removing the cut-off at 0.19 keV is that the DM contours extend further into the ${}^{8}\mathrm{B}$ region. However the signal region for these values of $m_{\mathrm{DM}}$ still lies between the ${}^{8}\mathrm{B}$ region and the ER band (where $pp$ neutrinos and other ER background events are expected to be detected). The lower right panel of Fig. 5 shows the median cross-section for LZ to make a discovery at $3\sigma$ (or greater) significance, while the inset shows a direct comparison of the ER efficiencies. For $m_{\mathrm{DM}}\sim 1~{}\mathrm{GeV}$, the different signal models result in similar sensitivity. The effect of removing the cut-off at 0.19 keV means that the efficiencies extend to small energies, while it allows the sensitivity of Models 2 and 3 to be enhanced with respect to Model 1 for $m_{\mathrm{DM}}\sim 0.1~{}\mathrm{GeV}$.

Appendix C Validation of the WS2013 background model

The profile likelihood ratio (PLR) test requires a background model. We here show that the simple ER background model that we have used gives a good fit to published LUX results. In particular, the background count rate and LUX’s background model for WS2013 were published in their axion search paper Akerib et al. (2017c). As in this work, only events measured with a radius smaller than 18 cm were considered. The measured data points and the LUX model are shown by the black data points and dotted black line in Fig. 6, as a function of S1 and $\log_{10}\mathrm{S2}$ in the upper and lower panels, respectively.

For our model of the WS2013 background, we assume a flat (in energy) component and a component from decays of ${}^{127}\mathrm{Xe}$, a cosmogenic isotope with a half-life of $\sim 36$ days that decayed during WS2013 (by WS2014-16, this component had completely decayed away). Decays of ${}^{127}\mathrm{Xe}$ lead to ER depositions at 5.2 keV, 1.1 keV and 0.19 keV with branching ratios of 13.1%, 2.9% and 0.6% respectively (we do not include the higher-energy decays that make up the remainder of the branching ratios).

The result of a fit of our model to the data is shown in Fig. 6. The yellow and blue lines show the contributions from the flat component and the ${}^{127}\mathrm{Xe}$ component, respectively. In both panels, we see that this simple model provides a good fit to the LUX data. It is also in good agreement with the published LUX model. The largest difference occurs in the $\log_{10}\mathrm{S2}$ comparison (lower panel) where we see that our background is slightly displaced to higher S2 values compared to the LUX model. However, our model is still in good agreement with the data points.

Appendix D Discovery potential with a higher background

In the main part of the paper, we assumed that LZ’s ER background will be dominated by pp and ${}^{7}\mathrm{Be}$ solar neutrinos. For a $5600\times 1000~{}\mathrm{kg\,day}$ exposure, we predict 252 events in the energy range from 1.5 to 6.5 keV, which is in good agreement with the value of 255 events predicted by the LZ Collaboration (LZ also include ${}^{13}\mathrm{N}$ solar neutrinos which could explain the small difference) Mount et al. (2017). In order for the solar neutrino background to dominate, the background contribution from dispersed radionuclides (particularly radon, krypton and argon) must be sub-dominant. This was the assumption made in Ref. Akerib et al. (2015). However, a more recent estimate in Ref. Mount et al. (2017) suggests that the rate from dispersed radionuclides could in fact dominate, resulting in a total of 1244 events for a $5600\times 1000~{}\mathrm{kg\,day}$ exposure in the energy range from 1.5 to 6.5 keV.

To assess the impact of a higher background on our LZ sensitivity projection, we recalculate the median cross-section for LZ to make a discovery at $3\sigma$ (or greater) significance assuming a background rate that is five times higher than the rate from only pp and ${}^{7}\mathrm{Be}$ solar neutrinos. This results in a total of 1260 events. The solar neutrino and radionuclide energy spectra are approximately flat in energy Schumann et al. (2015), so a simple rescaling of the pp and ${}^{7}\mathrm{Be}$ energy spectrum is a good approximation.

The resulting discovery cross-sections are shown in Fig. 7 for the more sensitive LZ scenario. The black solid line shows the result in the main part of the paper, where the ER background is dominated by pp and ${}^{7}\mathrm{Be}$ solar neutrinos, while the black dashed line shows the result when the ER background is five times larger. At low masses, the dark matter signal region is far from the ER backgrounds so the result does not change. At higher masses where the impact of the ER background is most significant, the discovery cross-section is only about a factor of two higher. This is because the PLR method still has some discrimination power between signal and background owing to the slight displacement of signal and background. Thus our projections for LZ’s sensitivity are robust against reasonable variations in the background rate.