Effective field theory search for high-energy nuclear recoils using the XENON100 dark matter detector

TL;DR

This paper re-analyzes XENON100 data using effective field theory to search for high-energy nuclear recoils from WIMPs, setting new limits on WIMP interactions and inelastic scattering models.

Contribution

It introduces a non-relativistic effective field theory framework to interpret XENON100 data for high-energy nuclear recoils, extending previous analyses.

Findings

No significant excess found over background.

Set 90% CL exclusion limits on WIMP-nucleon couplings.

Provided limits on inelastic WIMP scattering models.

Abstract

We report on WIMP search results in the XENON100 detector using a non-relativistic effective field theory approach. The data from science run II (34 kg $\times$ 224.6 live days) was re-analyzed, with an increased recoil energy interval compared to previous analyses, ranging from $(6.6 - 240)~\mathrm{keV_\mathrm{nr}}$. The data is found to be compatible with the background-only hypothesis. We present 90% confidence level exclusion limits on the coupling constants of WIMP-nucleon effective operators using a binned profile likelihood method. We also consider the case of inelastic WIMP scattering, where incident WIMPs may up-scatter to a higher mass state, and set exclusion limits on this model as well.

Tables

Table 1. Table 1: Definitions and contents of the analysis bins for the high energy channel. The expected background counts are calculated by taking the calibration sample and scaling it by 6.54 × 10 − 3 6.54 superscript 10 3 6.54\times 10^{-3} , which is the ratio of observed counts to calibration counts in a sideband.

#	Band	Energy Range $\left(\mathrm{cS1}\right)$	# Background Events	# Data Events
1	upper	30 - 40	24$\pm$5	20
2	upper	40 - 50	16$\pm$3	17
3	upper	50 - 80	12$\pm$3	11
4	upper	80 - 120	$1.1\pm 0.3$	1
5	upper	120 - 150	$\left(1.0\pm 0.5\right)\times {10}^{-1}$	1
6	upper	150 - 180	$\left(0.8\pm 0.4\right)\times {10}^{-1}$	0
7	lower	30 - 50	$0.9\pm 0.3$	0
8	lower	50 - 90	$\left(3.5\pm 1.2\right)\times {10}^{-1}$	0
9	lower	120 - 180	$\left(1.8\pm 0.7\right)\times {10}^{-1}$	0

Equations

$$\begin{array}{rl}& {\mathcal{O}}_1=1_{\chi }{1}_N\\ & {\mathcal{O}}_3=i{\vec{S}}_N\cdot (\frac{\vec{q}}{m_N}\times {\vec{v}}^{\perp })\\ & {\mathcal{O}}_4={\vec{S}}_{\chi }\cdot {\vec{S}}_N\\ & {\mathcal{O}}_5=i{\vec{S}}_{\chi }\cdot (\frac{\vec{q}}{m_N}\times {\vec{v}}^{\perp })\\ & {\mathcal{O}}_6=({\vec{S}}_{\chi }\cdot \frac{\vec{q}}{m_N})({\vec{S}}_N\cdot \frac{\vec{q}}{m_N})\\ & {\mathcal{O}}_7={\vec{S}}_N\cdot {\vec{v}}^{\perp }\\ & {\mathcal{O}}_8={\vec{S}}_{\chi }\cdot {\vec{v}}^{\perp }\end{array}\begin{array}{rl}& {\mathcal{O}}_9=i{\vec{S}}_{\chi }\cdot ({\vec{S}}_N\times \frac{\vec{q}}{m_N})\\ & {\mathcal{O}}_{10}=i{\vec{S}}_N\cdot (\frac{\vec{q}}{m_N})\\ & {\mathcal{O}}_{11}=i{\vec{S}}_{\chi }\cdot (\frac{\vec{q}}{m_N})\\ & {\mathcal{O}}_{12}={\vec{S}}_{\chi }\cdot ({\vec{S}}_N\times {\vec{v}}^{\perp })\\ & {\mathcal{O}}_{13}=i(\vec{S}\chi \cdot {\vec{v}}^{\perp })({\vec{S}}_N\cdot \frac{\vec{q}}{m_N})\\ & {\mathcal{O}}_{14}=i({\vec{S}}_{\chi }\cdot \frac{\vec{q}}{m_N})({\vec{S}}_N\cdot {\vec{v}}^{\perp })\end{array}$$

$${\mathcal{O}}_{15}=-({\vec{S}}_{\chi }\cdot \frac{\vec{q}}{m_N})\left[({\vec{S}}_N\times {\vec{v}}^{\perp })\cdot \frac{\vec{q}}{m_N}\right]$$

$$⟨cS1⟩=E_{\mathrm{nr}}\cdot (L_{\mathrm{y}}{\mathcal{L}}_{\mathrm{eff}})\cdot \left(\frac{S_{\mathrm{nr}}}{S_{\mathrm{ee}}}\right)$$

$$⟨{cS2}_{\mathrm{b}}⟩=E_{\mathrm{nr}}{\mathcal{Q}}_{\mathrm{y}}Y$$

$$p_{S2}({cS2}_{\mathrm{b}}|E)=\sum _{N^{\prime }}{P}_{\mathrm{pmt}}({cS2}_{\mathrm{b}}|Y{N}^{\prime },{\sigma }_Y\sqrt{N^{\prime }})\cdot \mathrm{Pois}(N^{\prime }|{\mu }_Q)$$

$$v_{\min }/c=\frac{1}{\sqrt{2m_N{E}_R}}\left|\frac{m_N{E}_R}{{\mu }_N}+{\delta }_m\right|$$

$${\sigma }_N^{\mathrm{SI}}={\left(C_1^N\right)}^2\frac{{\mu }_{\chi ,N}^2}{\pi }$$

$$\frac{\mathrm{d}R}{\mathrm{d}cS1}$$

$$=\int \frac{\mathrm{d}R}{\mathrm{d}E}G(cS1,E)\mathrm{d}E$$

$$R_{\mathrm{bi}{\mathrm{n}}_{\mathrm{i}}}=\frac{1}{2}{\int }_{{\mathrm{lower}}_i}^{{\mathrm{upper}}_i}\frac{\mathrm{d}R}{\mathrm{d}cS1}\mathrm{d}cS1$$

$$R_{\mathrm{bi}{\mathrm{n}}_{\mathrm{i}}}$$

$$=\int \frac{\mathrm{d}R}{\mathrm{d}E}{G}_i^{\prime }(E)\mathrm{d}E$$

Full text

Dark Matter, EFT, Xenon

I Introduction

Astrophysical and cosmological observations provide strong evidence that about 27% of the energy density of the universe is made out of Dark Matter (DM). The DM hypothesis is based on the existence of a non-luminous, non-baryonic, and non-relativistic particle, the nature of which is yet unknown Harvey et al. (2015); Bennett et al. (2013); Ade et al. (2016). Many well-motivated theoretical extensions of the Standard Model of particle physics predict the existence of one or more particles with the required properties, with masses and cross sections typically of the order of the weak scale. Such particles are collectively known as Weakly Interacting Massive Particles (WIMPs) Silk et al. (2010). The hypothesis that dark matter is constituted primarily of WIMPs is currently being tested by many experiments, either indirectly by seaching for evidence of their possible decay or annihilation in astrophysical processes, by searching for evidence of their direct production at collider experiments, or by directly measuring the rare scattering of astrophysical WIMPs from target nuclei in Earth-based laboratories Aprile et al. (2016); Tan et al. (2016); Akerib et al. (2016); Aalseth et al. (2013); Agnese et al. (2014); Angloher et al. (2012); Bernabei et al. (2010). We report on a search of this latter kind.

The traditional approach for computing predictions of the rate of WIMP-nucleon scattering has been to take only leading-order terms in a WIMP-nucleon effective field theory (EFT) with a very simple treatment of nuclear structure Lewin and Smith (1996). This leads to two main types of interactions, which are commonly labelled “Spin Independent” (SI) and “Spin Dependent” (SD). However, in recent years many authors have pointed out that in certain theories these interactions may be suppressed or nonexistent, such that otherwise subleading interactions may dominate the scattering process Chang et al. (2010). To account for this possibility in a systematic way, a more sophisticated EFT approach has been developed  Fitzpatrick et al. (2012); Anand et al. (2014); Fitzpatrick et al. (2013). In the new approach, an effective Lagrangian describing the WIMP-nucleus interaction is constructed, that takes into account all Galilean-invariant operators up to second order in the momentum exchange. This framework introduces new operators associated with different types of nuclear responses, along with the standard SI and SD ones, resulting in a set of fourteen operators $\mathcal{O}_{i}$ which may couple independently to protons and neutrons. In Eqs. (1) we list these operators following the convention from Anand et al. (2014). The operators depend explicitly on 4 linearly independent quantities: $\vec{v}^{\perp}\equiv\vec{v}+\frac{\vec{q}}{2\mu_{N}}$, the relative perpendicular velocity between the WIMP and the nucleon, $\vec{q}$, the momentum transferred in the scattering event, and $\vec{S}_{\chi}$, $\vec{S}_{N}$, the WIMP and nucleon spins. $\mathcal{O}_{2}$ is not considered here as it cannot be obtained from a relativistic operator at leading order.

[TABLE]

[TABLE]

Unlike the more commonly studied types of interaction (SI,SD), which are not suppressed when $\vec{q}\rightarrow 0$ and for which the scattering rate on nucleons is expected to be largest for low energy nuclear recoils, some of the new EFT operators depend explicitly on $\vec{q}$ and so their interaction cross section is suppressed for low momentum transfers. Consequently, their scattering rate peaks at non-zero nuclear recoil energy. For sufficiently high WIMP masses, this may even occur outside typical analysis windows, which usually have an upper range of around $43\,\mathrm{keV_{\mathrm{nr}}}$ (nuclear recoil equivalent energy) since they are designed to search for SI and SD interactions, which predict exponentially-falling recoil spectra (see Figure 1). Due to the theoretical bias of only considering SI and SD interactions, high energy nuclear recoils remain unexplored in many experiments.

Another typical assumption that can be relaxed is that WIMPs should scatter elastically with nuclei. There exist dark matter models in which the incoming and outgoing WIMPs have different mass states Tucker-Smith and Weiner (2001) separated by a keV-scale splitting. In the case where the outgoing state is more massive than the incoming state, the cross section for low recoil energies can again be suppressed, this time by scattering kinematics. Recently an inelastic adaptation of the EFT operator framework discussed above was developed Barello et al. (2014). In this case the operators presented in Eqs. 1 are modified such that $\vec{v}^{\perp}_{inelastic}=\vec{v}^{\perp}_{elastic}+\frac{\delta_{m}}{|{\vec{q}}|^{2}}\vec{q}$. We consider this case in section III.3.2.

The EFT framework of Fitzpatrick et al. (2012) is constructed at the WIMP-nucleon level and so each operator may be present independently for protons and neutrons, though UV models can of course correlate their couplings. The full EFT thus has 28 coupling parameters in addition to the WIMP mass, plus a mass splitting $\delta$ in the inelastic case. This parameter space is too large to explore in full, so we take a similar approach to the SI/SD case

and assume only one active operator at a time, considering it equally coupled to protons and neutrons (the “isoscalar” case). However, to facilitate the full exploitation of these results by the community, we provide in supplementary material a set of tools for converting any theoretical recoil spectrum $\mathrm{d}R/\mathrm{d}E$ into an accurate event rate prediction for this analysis, including all detector response and analysis efficiency effects. This may help to set a mildly conservative but quite accurate limit on arbitrary models in the full EFT parameter space, or any other particle dark matter model for which one can supply the expected recoil spectrum. These tools are described further in Appendix A.

Motivated by these EFT extensions of the standard WIMP framework, we report on an analysis extending the searched recoil energy range up to $240~{}\mathrm{keV_{\mathrm{nr}}}$ for the first time in the XENON100 experiment, and present exclusion limits on all operators for both elastic and inelastic WIMP cases.

II The XENON100  Detector

The XENON100 detector is a cylindrical dual-phase xenon (liquid and gas) time projection chamber (TPC). It is installed at the Laboratori Nazionali del Gran Sasso (LNGS) in Italy and contains 161 kg of liquid xenon (LXe), of which 62 kg function as the active target  Aprile et al. (2012a). The detector uses of a total of 178 1-inch square Hamamatsu R8520-AL photomultiplier tubes (PMTs) employed in two arrays, one in the gas phase at the top of the TPC, and the other at the bottom, immersed in the LXe.

A particle interacting with the LXe deposits energy that creates both prompt scintillation ($\mathrm{S1}$) and delayed proportional scintillation ($\mathrm{S2}$) which are detected using the two PMT arrays. The $\mathrm{S2}$ signal is produced by ionization electrons, drifted in an electric field of $530$V/cm towards the liquid-gas interface, where they are extracted to the gas phase using a stronger electric field of $\sim 12$kV/cm in which the proportional scintillation occurs. The spatial distribution of the $\mathrm{S2}$ signal on the top PMT array, together with the time difference between $\mathrm{S1}$ and $\mathrm{S2}$ signals, provide respectively $x$-$y$ and $z$ position information for each interaction, allowing 3D position reconstruction to be achieved.

Interaction in different locations of the detector have different signatures. In order to take these effects into account, a correction is applied based on light and charge collection efficiency maps. These maps are prepared using calibration sources ranging up to energies well above $240~{}\mathrm{keV_{\mathrm{nr}}}$, which is the highest energy recoil considered in this paper. The corrected signals ($\mathrm{cS1}$,$\mathrm{cS2}_{\mathrm{b}}$) are spatially independent and uniform to all interactions Aprile et al. (2012a). Note that some of the top PMTs saturate for large $\mathrm{S2}$ signals and we therefore use in this analysis only the bottom PMT array to infer the energy scale in $\mathrm{S2}$.

The $\mathrm{S1}{}/\mathrm{S2}{}$ ratio is known to differ between nuclear recoil (NR) and electronic recoil (ER) interactions, and is thus used as a discriminating variable between a WIMP signal and ER background. The logarithm of this ratio, $\log(\mathrm{cS2}_{\mathrm{b}}{}/\mathrm{cS1}{})$ is referred later in the text as the discriminating “$y$” variable.

III Data Analysis

In this work we re-analyze science run II data recorded between February 2011 and March 2012, corresponding to 224.6 live days. The characterization of the detector response to ER interactions is performed using dedicated calibration campaigns with 60Co and 232Th radioactive sources, while the response to NR interactions is performed using 241AmBe neutron source calibration campaigns.

This work extends the previous results Aprile et al. (2012b, 2016), referred to in the following as the low-energy channel, with a new study exploring the recoil energy range between $43-240~{}\mathrm{keV_{\mathrm{nr}}}$. The data analysis is divided into two mutually exclusive channels, one optimized for low energies and ranging from 3-30 PE in $\mathrm{cS1}$ (low-energy), the other optimized for high energies recoils ranging from 30-180 PE in $\mathrm{cS1}$ (high-energy). These two analyses are then combined statistically.

III.1 Low energy channel

This analysis channel relies on the re-analysis of run II data described in Aprile et al. (2016). The region of interest (ROI), the background expectation models, data selections and their acceptances are mostly unchanged and so are only briefly summarized here. Differences with respect to said results are highlighted when present.

The ROI for this channel is defined in the ($y,\mathrm{cS1}{}$)-plane and is shown in Figure 2. The lower bound on $y$ corresponds to a 3 $\sigma$ acceptance quantile (as a function of $\mathrm{cS1}$) of a 20 GeV WIMP mass signal model assuming an $\mathcal{O}_{1}$ (SI) interaction, while the upper bound is fixed at $y=2.7$. The range in $\mathrm{cS1}$ is selected as 3 to 30 PE. The ROI is further divided into eight sub-regions (also called bands) depending on the operator $\mathcal{O}_{i}$ and on the WIMP mass hypothesis. These bands are arranged to achieve constant expected signal density in each region, as described in Aprile et al. (2016).

Other than falling into the ROI, an event should fulfill several additional selection criteria (cuts). Data quality and selection cuts are defined to remove events with poor data quality or noisy signals. Events are discarded if they present a time-coincident signal in the outer LXe veto, $\mathrm{S2}$ signals below threshold, multiple-scatters, or are localized outside a predefined fiducial volume of 34 kg. In addition, this analysis channel uses the post-unblinding cuts and data reprocessing described in Aprile et al. (2016). More details on these selection criteria and their relative WIMP signals acceptances can be found in Aprile et al. (2014a, 2016).

Note that this analysis channel does not employ a variable lower $\mathrm{S1}$ threshold as a function of the event position in the TPC, but instead applies a fixed lower threshold cut on $\mathrm{cS1}$ at 3 PE, conversely to the choice made in Aprile et al. (2016).

The expected background is modeled separately for ER and NR contributions which are then scaled to exposure and added together. The NR background is estimated by Monte Carlo simulation and accounts for the radiogenic and cosmogenic neutron contributions Aprile et al. (2013a). The ER background is parametrized as the linear combination of Gaussian-shaped and non-Gaussian components. The former is obtained via a parametric fit of the 60Co and 232Th calibration data, as discussed in Aprile et al. (2012b).

The latter, which consist of anomalous events such as those presenting incomplete charge collection or accidental coincidence of uncorrelated $\mathrm{S1}$s and $\mathrm{S2}$s, is evaluated via dedicated techniques described in Aprile et al. (2016).

Systematic uncertainties on the background model arising from the Gaussian parametrized fit, and from the normalisations of the NR and non-Gaussian components, have been evaluated and propagated to each band. These errors are small with respect to the statistical uncertainties of each band, which are conservatively taken as the overall uncertainty Aprile et al. (2016), as discussed in Sec. III.4.

III.2 High energy channel

This analysis channel targets high energy nuclear recoils and is the focus of this work. The data selection criteria used are based on the criteria described in detail in Aprile et al. (2014a), which were optimized for high acceptance to low energy nuclear recoils. Most of these cuts were found to be fully compatible with (or easily extended) to high energy depositions, however some required more comprehensive studies, which are described in the following .

The width of an $\mathrm{S2}$ pulse increases with the depth (z) of the interaction. This is due to the diffusion of the electron cloud during its propagation through the liquid xenon. Since low energy $\mathrm{S2}$ events show larger spread due to low statistics of drifted electrons, the cut was previously defined in an energy-dependent way. However, for the large recoil energies considered in this channel, this energy dependency is no longer valid. We therefore use here a cut on the $\mathrm{S2}$ width which is a function of the depth of the interaction alone.

As a WIMP will interact only once in the detector, we remove events which have more than one $\mathrm{S2}$. We adopt in this analysis a cut that is more suitable to higher energies and demand a single $\mathrm{S2}$ in a 160 $\mu$s window, instead of a linear dependence between the second $\mathrm{S2}$ size and the first.

To define the interaction’s exact location in ($x,y$), we use several algorithms, one of which is based on a Neural Network (NN) Aprile et al. (2014a). The NN was not trained to recognize high energy ER events and therefore a cut on the NN reconstruction quality is not suitable for this analysis. We therefore discard this cut but keep all other selections on position reconstruction quality, which is sufficient to ensure a correct position reconstruction.

The total acceptance to WIMP signals is computed based on 241AmBe calibration data as a function of $\mathrm{cS1}$, following the procedure described in Aprile et al. (2014a). We present this function in Figure 3, where the total acceptance is fitted using a third order polynomial.

We define our signal region in the discrimination $(y,\mathrm{cS1})$-plane using 241AmBe calibration data. The region of interest is shown in Figure 2 as blue contour lines. The upper bound in $y$ is defined such that the contribution due to xenon inelastic interaction lines is negligible. The lower bound is defined as the 3 $\sigma$ acceptance quantile of the 241AmBe distribution.

We divide our signal region into two bands in $y$, constructed such that the 241AmBe data sample is equally distributed in between them. The number of events in each band is $\sim 3000$. The bands are further divided into nine bins, the number and boundaries of which have been optimized via Monte-Carlo (MC) simulation. The definitions of the bins boundaries are presented in Table 1 and in Figure 2.

The main source of background results from ER leakage. We therefore estimate the background distribution in the ROI using 60Co and 232Th calibration events. Contributions from radiogenic and cosmogenic neutrons, as well as accidental coincidence, are negligible for such a high energy recoil. In Table 1 we report the background expectation in the ROI along with the observed events for each bin. Here the background expectation is computed by scaling the calibration sample yield by $6.54\times 10^{-3}$, which is the ratio of observed counts to calibration counts in an independent sideband. The sideband is defined above the upper limit of this analysis and below the ER calibration band mean. Note that in the computation of exclusion limits, the background normalization is fitted to data, rather than using the sideband normalization, as described in section III.4.

III.3 Signal model

The signal model is produced by taking a theoretical event rate spectrum, the production of which is described in sections III.3.1 and III.3.2, and applying the analysis acceptance and detector response as described in  Aprile et al. (2014a) to obtain the expected event rate in the detector in terms of detector variables (i.e. $\mathrm{cS1}$, $\mathrm{cS2}_{\mathrm{b}}$). In both analysis channels, we use Eq. 2 in order to compute the expected average $\mathrm{cS1}$ for a given NR energy,

[TABLE]

where $E_{\mathrm{nr}}$ is the recoil energy, $L_{\mathrm{y}}$ is the average light yield in the detector, $\mathcal{L}_{\mathrm{eff}}$ is the scintillation efficiency relative to 122$\mathrm{keV_{\mathrm{ee}}}$ as a function of $E_{\mathrm{nr}}$, and $S_{\mathrm{ee}}$ and $S_{\mathrm{nr}}$ are the quenching factors due to the externally applied electric field. Aside from $E_{\mathrm{nr}}$ and $\mathcal{L}_{\mathrm{eff}}$ these parameters have fixed values, namely $L_{\mathrm{y}}=2.28\pm 0.04$, $S_{\mathrm{nr}}=0.95$, and $S_{\mathrm{ee}}=0.58$. Recoils below $3~{}\mathrm{keV_{\mathrm{nr}}}$ are assumed to produce no light. For details of the physics behind these parameters and the construction of the signal probability density function (PDF) please see Aprile et al. (2014a, 2016).

For the low-energy region, the expected $\mathrm{cS2}_{\mathrm{b}}$ signal is computed following Aprile et al. (2013b) using Eq. 3,

[TABLE]

where $Y=8.3\pm 0.3$ is the amplification factor determined from the detector response to single electrons Aprile et al. (2014b), and $\mathcal{Q}_{\mathrm{y}}$ is the charge yield as a function of $E_{\mathrm{nr}}$. Applying the detector and PMT responses, and the acceptance as in Aprile et al. (2016), defines the low-energy signal model over the region $3~{}\mathrm{PE}<\mathrm{cS1}{}<30~{}\mathrm{PE}$, with $\mathrm{cS2}_{\mathrm{b}}{}>73.5~{}\mathrm{PE}$ as the $\mathrm{S2}$ threshold.

Eq. 3 hides a subtlety. The actual $\mathrm{cS2}_{\mathrm{b}}$ PDF is composed of two pieces, a Poisson term associated with the initial charge liberation and a Gaussian term associated with the PMT response and other detector effects:

[TABLE]

where $\mu_{Q}=E_{\mathrm{nr}}\mathcal{Q}_{\mathrm{y}}$ is the expected number of liberated charges in a nuclear recoil event of energy $E$, and $N^{\prime}$ is the actual number of liberated charges. The amplification factor $Y$ is applied to the actual number of liberated charges $N^{\prime}$, not the expected number $\mu_{Q}$. Associated with this is the variance of the Gaussian response PDF, $\sigma_{Y}\sqrt{N^{\prime}}$, where in this analysis $\sigma_{Y}=6.93$ as measured and described in Aprile et al. (2014b).

For the high energy region we cannot produce the $\mathrm{S2}$ distribution in the same way as the method in Aprile et al. (2013b), since it has not been calibrated for such high recoil energies. We therefore use the NR calibration data distribution in log($\mathrm{\mathrm{cS2}_{\mathrm{b}}/\mathrm{cS1}}$) to estimate the WIMP distribution. Above 180 PE in $\mathrm{cS1}$, the event yield of 241AmBe data is too low to estimate the distribution accurately. This forms the upper bound of this analysis. With the $\mathrm{cS2}_{\mathrm{b}}$ distribution determined by this empirical method, we require only a prediction of the $\mathrm{cS1}$ distribution. This is obtained from Equation (2), followed by the application of detector and PMT responses, as well as the acceptance given in Figure 3, which completes the high-energy signal model definition.

Figures  4 and 5 shows dashed signal distribution examples for two EFT operators and for the low and the high energy region, respectively. In both cases, the signal distributions are normalized to yield 5 events in the total energy range (low-energy and high-energy).

III.3.1 Elastic scattering

The expected recoil energy spectrum of each WIMP mass for each EFT operator is calculated using the Mathematica package DMFormFactor supplied by Anand et. al. Fitzpatrick et al. (2013); Anand et al. (2014). We use standard assumptions as in previous analyses (e.g Aprile et al. (2016)) regarding the local dark matter density and velocity distribution, namely $\rho_{\mathrm{local}}=0.3$ GeV$\cdot c^{-2}$/$\mathrm{cm}^{3}$ and a Maxwell-Boltzman distribution with a mean given by the local circular velocity $v_{0}=220$ km/s and cut off at an escape velocity of $v_{\mathrm{esc}}=544$ km/s. The responses of xenon nuclei to a scattering event are computed from one-body density matrices provided with the package, in contrast to the Helm form factors which have been used in previous analyses. These spectra are produced for the seven most abundant xenon isotopes (128, 129, 130, 131, 132, 134 and 136), combined in proportion to the abundance of these isotopes in the XENON detector Aprile et al. (2013c), then translated into expected signal rates via the method described above.

III.3.2 Inelastic WIMP scattering

To obtain recoil spectra for WIMP-nucleon scattering for all EFT operators with inelastic kinematics, we use a modified version of DMFormFactor provided by Barello et. al. Barello et al. (2014). The authors have modified the original package to enforce the new energy conservation condition $\delta_{m}+\vec{v}\cdot\vec{q}+\left|\vec{q}\right|^{2}/2\mu_{N}=0$, primarily by replacing $\vec{v}^{\perp}_{elastic}\rightarrow\vec{v}^{\perp}_{inelastic}=\vec{v}^{\perp}_{elastic}+\frac{\delta_{m}}{|{\vec{q}}|^{2}}\vec{q}$ in the definitions of the EFT and nuclear operators, giving rise to the well-known minimum velocity for scattering

[TABLE]

where $\mu_{N}$ is the WIMP-nucleon reduced mass.

Assumptions regarding the dark matter halo and nuclear physics are unchanged. The mass splitting $\delta_{m}$ between dark matter states is varied from 0 to 300 keV, safely beyond the value at which the predicted rate is zero for the entire mass range we consider.

III.4 Statistical inference

The statistical interpretation of data is performed using a binned profile likelihood method, in which hypothesis testing relies upon a likelihood ratio test statistic, $\tilde{q}$, and its asymptotic distributions Cowan et al. (2011). The two analysis channels are combined by multiplying their likelihoods together to produce a joint likelihood. Both analyses parametrize the NR relative scintillation efficiency, $\mathcal{L}_{\mathrm{eff}}$, based on existing measurements Aprile et al. (2011a). Its uncertainty is the major contributor to energy scale uncertainties and is considered as correlated between the two analysis channels via a joint nuisance likelihood term. Throughout this study, all the parameters related to systematic uncertainties are assumed to be normally distributed.

For the low energy channel an extended likelihood function is employed which is very similar to the one reported in Aprile et al. (2011b) and described in detail in Aprile et al. (2016). The $(y,\mathrm{cS1}{})$-plane is divided into eight WIMP mass dependent bands where events are counted. This binned approach is extended with the corresponding $\mathrm{cS1}$-projected PDF of each band. The total normalization of the background is fit to data, and an uncertainty is assigned to the relative normalization of each band according to the corresponding statistical uncertainty of the calibration sample. Signal shape variations due to energy scale uncertainty are modeled via simulation. These include the said $\mathcal{L}_{\mathrm{eff}}$ uncertainties and additionally the charge yield uncertainties, which are parametrized based on $\mathcal{Q}_{\mathrm{y}}$ measurement as described in Aprile et al. (2013b).

The high energy channel analysis employs a binned likelihood function. Observed and expected event yield are compared in the nine ROI $(y,\mathrm{cS1}{})$-bins described in section III.2. Given the large statistical uncertainty of the background model the above extended likelihood approach is not repeated here. Instead, the maximum likelihood estimation of the background expectation in each bin is constrained by the statistical uncertainty of the calibration sample, while the total normalization is fit to the data. Additionally, to account for potential mismodeling of the expected background distribution, mainly due to anomalous multiple scatter events, a systematic uncertainty of 20% is assigned independently to each bin. In the high energy channel, uncertainty on the signal acceptance of analysis selections are computed for each signal hypothesis using the parametrized acceptance curve shown in Figure 3. Uncertainties on the signal model $(y,\mathrm{cS1}{})$ distribution due to 241AmBe sample statistical fluctuations, as well as energy scale shape variation due to $\mathcal{L}_{\mathrm{eff}}$ uncertainties, are taken into account.

IV Results

A benchmark region of interest is defined between the upper and lower thresholds in $\mathrm{cS1}$ for each channel. This region is bounded in $y$-space from above by the 241AmBe NR mean line and below by the lower 3$\sigma$ quantile of the 241AmBe neutron calibration data. The expected background in the region is $3.0\pm 0.5_{stat}$ (low-energy) and $1.4\pm 0.3_{stat}$ (high-energy). The number of DM candidates in this benchmark region is 3 (low-energy), and 0 (high-energy). Consequently, the data is compatible with the background-only hypothesis and no excess is found.

For the elastic scattering case, a 90% CLS Read (2000) confidence level limit is set on the effective coupling constant, $c_{i}$, for all operators and masses in the range of 10 GeV/$c^{2}$ to 1 TeV/$c^{2}$. The $c_{i}$ are dimensionful, with units of $[\mathrm{mass}]^{-2}$, so we first convert them to dimensionless quantities by multiplying them by $m_{\mathrm{weak}}^{2}=(246.2\text{ GeV})^{2}$, following the conventions of Anand et al. (2014).

These limits are shown in Fig. 6 in black, along with limits from CDMS-II Si, CDMS-II Ge and SuperCDMS Schneck et al. (2015).

For the inelastic scattering case, 90% CLS confidence level limits on the coupling constants (again scaled by $m_{\mathrm{weak}}^{2}$) are set. Fig. 7 shows limits on the $\mathcal{O}_{1}$ (SI) coupling constant as a function of mass splitting and WIMP mass, Fig. 8 shows limits for all other operators as a function of the mass splitting $\delta_{m}$ with a fixed WIMP mass of 1 TeV/$c^{2}$, projections of results from CDMS-II Ahmed et al. (2010), ZEPLIN-III Akimov et al. (2010), and XENON100 Aprile et al. (2011c) in the coupling constant and $\delta_{m}$ parameter space are also reported.

For the elastic operator $O_{1}$ our results can be compared to those of standard SI analyses by computing the relevant zero-momentum WIMP-nucleon cross-sections. This is not simple to do rigorously because the treatment of nuclear structure used in our analysis is different than in standard analyses, however this difference is small for scattering via $O_{1}$. We can therefore quite safely use the ‘traditional’ correspondence De Simone and Jacques (2016)

[TABLE]

where $\mu_{\chi,N}$ is the WIMP-nucleon reduced mass. Standard SI analyses assume isospin-conserving interactions, as we do in this analysis, so we can simply set $C^{N}_{1}=C^{0}_{1}$, such that $\sigma_{p}^{\mathrm{SI}}=\sigma_{n}^{\mathrm{SI}}$.

In principle a similar comparison can be done between our limit on the $O_{4}$ coupling and standard SD analysis limits, however this time the standard analyses do not assume isospin-conserving interactions. Instead they typically assume maximal isospin violation, that is, assuming that WIMPs couple either protons or neutrons. Limits are then derived independently on $\sigma_{p}^{\mathrm{SD}}$ and $\sigma_{n}^{\mathrm{SD}}$. Because of this difference in assumptions, our limits on SD couplings are not directly comparable to usual analyses. However, they can be recast under the appropriate alternate model assumptions using the detector response tables we provide in the supplementary material.

V Summary

We have shown the first analysis of XENON100 data at recoil energies above 43 keV, with the new high energy bound set to 240 keV. We considered in this paper two models which predict interactions in this energy region: an EFT approach for elastic WIMP-nucleon scattering, and a similar EFT approach but considering instead inelastic WIMP-nucleon scattering. The observed data was compatible with background expectations, and 90% CLS exclusion limits were constructed for WIMP masses between 10-1000 GeV.

Acknowledgements.

We would like to thank Andrew Liam Fitzpatrick and Spencer Chang for supplying and helping with their Mathematica packages . We gratefully acknowledge support from the National Science Foundation, Swiss National Science Foundation, Deutsche Forschungsgemeinschaft, Max Planck Gesellschaft, German Ministry for Education and Research, Netherlands Organisation for Scientific Research, Weizmann Institute of Science, I-CORE, Initial Training Network Invisibles (Marie Curie Actions, PITNGA-2011-289442), Fundacao para a Ciencia e a Tecnologia, Region des Pays de la Loire, Knut and Alice Wallenberg Foundation, Kavli Foundation, and Istituto Nazionale di Fisica Nucleare. We are grateful to Laboratori Nazionali del Gran Sasso for hosting and supporting the XENON project.

Appendix A SIGNAL MODEL DETECTOR RESPONSE TABLE

In this appendix we describe digital tables which can be used to construct an accurate signal model for this analysis given any input recoil spectrum $\mathrm{d}R/\mathrm{d}E$ arising from a theoretical model. A visualization of the tables is shown in Fig. 9, and in section A.1 we show a simple example Python code of how to use the supplied tables. Currently we provide these tables only for the high-energy analysis region.

The signal model for the high-energy analysis region can be expressed analytically in the form:

[TABLE]

where $\epsilon_{\mathrm{S1}}(\mathrm{cS1})$ and $\epsilon_{\mathrm{S2^{\prime}}}(E)$ represent analysis cut efficiencies, $p_{\mathrm{S1}}(\mathrm{\mathrm{cS1}}|E)$ encodes detector effects, and $\mathrm{d}R/\mathrm{d}E$ gives the theoretically predicted nuclear recoil rate from WIMP scattering. In the second line we emphasis that all the detector and analysis effects can be encoded in a single function $G(\mathrm{cS1},E)$. To make a signal prediction for the bins in our analysis, this expression needs to be integrated over the appropriate range of $\mathrm{cS1}$ for each bin (and divided by two to account for the banding structure in $\mathrm{cS2}_{\mathrm{b}}$):

[TABLE]

With some simple rearrangement this rate can be written in terms of an integral over the detector response function $G$ as follows

[TABLE]

where in the last line we absorb the factor of $1/2$ into the definition of $G^{\prime}_{i}$. We see here that the signal rate for each bin can be expressed as an integral over the recoil spectrum times a detector response function $G^{\prime}_{i}$ for that bin. It is these detector response functions which are shown in Fig. 9, and which we provide digitally for use by the community. A low-resolution example is given in Table 2. With these tables it is simple to produce a signal model for our analysis for any theoretical recoil spectrum. The functions $G^{\prime}_{i}$ are provided for three values of the nuisance variable $\mathcal{L}_{\mathrm{eff}}$, namely the median value and values at $\pm 1\sigma$ in $\mathcal{L}_{\mathrm{eff}}$. From these, along with the measured background rates given in table 1, one may construct a likelihood which accounts for uncertainties in $\mathcal{L}_{\mathrm{eff}}$, Alternatively simply using the $-1\sigma$ value produces quite an accurate prediction and is generally conservative.

A.1 Example code

⬇

import numpy as np

from numpy import newaxis

from scipy.interpolate import interp1d

def TrapI(x,y):

"""Simple trapezoid integration"""

w = x[1:] - x[:-1]

h = (y[1:] + y[:-1])/2.

return np.sum(w*h,axis=0)

Load detector response table

data = np.loadtxt("detector_table.dat")

E = data[:,0]; Gi = data[:,1:]

Load test recoil spectrum (1 TeV WIMP, O6)

data = np.loadtxt("O6_1TeV.dat")

Er = data[:,0]

Input spectra is normalised to coupling^2=1,

rescale to something near limit (1e3)

Also multiply in the appropriate exposure

dRdE = data[:,1] * (1e3/1.) * 224.6*34.

Interpolate recoil spectrum to table values

Assume spectrum zero outside data given

f_dRdE = interp1d(Er,dRdE)

dRdE_matched = f_dRdE(E)

Ri = TrapI(E[:,newaxis],Gi*dRdE_matched[:,newaxis])

for i,R in enumerate(Ri):

print "bin {0}: rate = {1:.2g}".format(i+1,R)

Output:

bin 1: rate = 0.081

bin 2: rate = 0.098

bin 3: rate = 0.35

bin 4: rate = 0.46

bin 5: rate = 0.29

bin 6: rate = 0.22

bin 7: rate = 0.18

bin 8: rate = 0.47

bin 9: rate = 0.84