Prospects for dark matter detection with inelastic transitions of xenon

Dark matter can scatter and excite a nucleus to a low-lying excitation in a direct detection experiment. This signature is distinct from the canonical elastic scattering signal because the inelastic signal also contains the energy deposited from the subsequent prompt de-excitation of the nucleus. A measurement of the elastic and inelastic signal will allow a single experiment to distinguish between a spin-independent and spin-dependent interaction. For the first time, we characterise the inelastic signal for two-phase xenon detectors in which dark matter inelastically scatters off the Xe-129 or Xe-131 isotope. We do this by implementing a realistic simulation of a typical tonne-scale two-phase xenon detector and by carefully estimating the relevant background signals. With our detector simulation, we explore whether the inelastic signal from the axial-vector interaction is detectable with upcoming tonne-scale detectors. We find that two-phase detectors allow for some discrimination between signal and background so that it is possible to detect dark matter that inelastically scatters off either the Xe-129 or Xe-131 isotope for dark matter particles that are heavier than approximately 100 GeV. If, after two years of data, the XENON1T search for elastic scattering nuclei finds no evidence for dark matter, the possibility of ever detecting an inelastic signal from the axial-vector interaction will be almost entirely excluded.


Introduction
The properties of particle dark matter remain unknown. Searches with direct detection experiments are one of the most promising ways of detecting dark matter through an interaction other than gravity. A positive detection is expected to yield information on the particle mass, the cross-section and information on the form of the interaction [1,2]. Although there has not yet been a conclusive detection [3][4][5][6][7][8][9][10], direct detection experiments have demonstrated a remarkable record of increasing their sensitivity by an order of magnitude approximately every three years and this increase is expected to continue over the next decade [11]. Two-phase xenon experiments have proven to be particularly sensitive and we are approaching the tonne-scale era with the LUX [12] and XENON1T [13] experiments, and funding has been secured for the approximately five-tonne successor experiments, LZ [14,15] and XENONnT [16,17]. There is also a longer-term proposal for DARWIN [18], an even larger ∼ 20-tonne experiment whose aim is to explore all of the dark matter parameter space not limited by neutrino backgrounds [19]. The left diagram shows the canonical elastic scattering process where the dark matter simply causes the nucleus to recoil; an experiment measures the number of events and the nuclear recoil energy. The right diagram depicts the inelastic scattering process. In this case, the dark matter excites the xenon isotope which then promptly decays emitting a photon. For the 129 Xe and 131 Xe isotopes of xenon, the photon/excitation energies are 39.6 keV and 80.2 keV respectively. We assume that the photon mean free path is sufficiently short that the experiment measures the recoil of the nucleus at the same time as the prompt de-excitation photon.
Multi-tonne xenon experiments bring new opportunities to search for rare signals. This is for two reasons. Firstly, the larger target mass means that there are more xenon nuclei for the dark matter to scatter off and secondly, larger experiments allow for backgrounds to be significantly reduced, even down to the irreducible background from coherent neutrino scattering. This is because more of the liquid xenon can be used to self-shield the fiducial volume where dark matter signals are searched for.
The canonical search with direct detection experiments is the elastic scattering process depicted in the left diagram of figure 1. The interaction with the dark matter particle causes the xenon nucleus to recoil with an energy typically in the range 1-100 keV. Since some nuclear isotopes have excitations in this energy range, it was long ago realised that these nuclear excitations could also play a role in the detection of dark matter [20,21]. In this case, some part of the energy transferred from the dark matter particle causes the excitation of the nucleus while the other part causes the nucleus to recoil. The excited nucleus then decays emitting a photon. This process is depicted in the right panel of figure 1. For experiments with xenon, there are two isotopes of interest, 129 Xe and 131 Xe, which make up 26.4% and 21.2% of natural xenon and have an excitation energy and lifetime of 39.6 keV and 80.2 keV, and 0.97 ns and 0.48 ns respectively. In this process the recoil energy of the nucleus and the energy from the prompt de-excitation of the nuclear isotope are measured. The experimental resolution of xenon detectors is ∼ 10 ns [22] so the short lifetimes mean that the recoil and de-excitation cannot be separately resolved. However, the mean free path of the de-excitation photon is O(1) mm [23], comparable to the spatial resolution of xenon detectors [22], so a dedicated pulse-shape analysis may partially resolve the nuclear recoil and the photon energy deposition for a fraction of the events. We leave an analysis of the pulse-shape for the future and here make the assumption that the detector cannot resolve the recoil energy and photon energy i.e. it is only the total energy that is measured.
Nuclear structure functions are required in order to accurately predict the cross-section for dark matter to excite the nucleus. It is only recently that precision shell-model calculations of the structure functions for xenon isotopes have become available [24][25][26] (see also [27,28] for earlier calculations). The contribution of different nucleons to inelastic scattering does -2 -

JCAP05(2016)033
not add coherently so this process does not benefit from the nucleon-number-squared (∼ 10 4 ) enhancement of elastic spin-independent interactions [26]. The current absence of any elastic signal means that for these interactions, experiments would have to improve their sensitivity by at least this factor to begin to see the inelastic signal. Such a large sensitivity gain is unlikely to be achieved in the foreseeable future. In contrast, the structure functions for elastic and inelastic processes are more comparable for the axial-vector interaction, with the elastic structure function being only around 10 times larger [25]. This is because elastic scattering for the axial-vector interaction is spin-dependent and also does not have the the nucleon-numbersquared enhancement [24]. The initial discovery of dark matter will not be made with the inelastic process for the axial-vector interaction (because of the additional suppression of the inelastic rate from the structure function and the scattering kinematics). However, a detection of the inelastic signal would provide additional information to complement the elastic scattering signal. As a trivial example, while the elastic signal may come from either a spin-independent or spin-dependent interaction, the inelastic signal will only be detectable for a spin-dependent interaction so it detection would strongly point to a spin-dependent interaction. Further implications of measuring the inelastic signal are left for a future paper.
A number of single-phase xenon experiments have searched for the 39.6 keV de-excitation from the 129 Xe isotope [29][30][31] (see also [32]). However these experiments generally set weaker limits than two-phase xenon experiments because they do not have the same ability as twophase experiments to discriminate between signal and background processes. No search or sensitivity study has been carried out for a two-phase xenon detector. This is the aim of this paper: to characterise the inelastic scattering signal for a two-phase xenon detector, quantify the sensitivity of upcoming tonne-scale experiments to this inelastic process and assess whether a future detection can be made.
Our paper is structured as follows. In section 2 we recap the basic principles of dark matter scattering and describe how we model the elastic and inelastic signals in two benchmark xenon detectors. In section 3 we discuss the main backgrounds and calculate both the signal and background distributions in terms of the parameters that a xenon detector measures. Section 4 describes our frequentist method for calculating the sensitivity of upcoming tonne-scale experiments while section 5 contains our main results. We end with a discussion of interesting follow-up studies and our conclusions in section 6. A number of short appendices gather the formulae that we use for the generation of photons and electrons for nuclear and electronic interactions, a check of the statistical method that we employ, an explicit demonstration that the LUX neutron-only limits are generally stronger than the PICO proton-only limits and finally, a check of our results under an alternative dark matter halo model.

Modelling elastic and inelastic recoils of xenon
In this section we first review the usual formalism for elastic scattering of dark matter with a xenon nucleus in terms of the recoil energy of the nucleus. We show that this is easily extended to the case of inelastic scattering. Xenon detectors do not directly measure the energy but rather the scintillation light. We describe our modelling of the generation and detection of the scintillation light, which is based on the NEST formalism [33][34][35][36]. We then describe the properties of present and upcoming tonne-scale direct detection experiments and discuss the observable signals and their rate.

Scattering rates
The differential event rate for both elastic and inelastic scattering of dark matter with a xenon nucleus of mass m A in the detector frame may be written as where E R is the recoil energy of the xenon nucleus, m DM is the dark matter mass, ρ DM = 0.3 GeV/cm 3 is the local dark matter density [37], v and v are the dark matter speed and velocity, and f DM ( v) is the dark matter velocity distribution in the galactic frame. We assume the isothermal Standard Halo Model so that f DM ( v) ∝ exp −v 2 /v 2 0 is a Maxwell-Boltzmann distribution in the galactic frame with a hard cut-off at the galactic escape speed v esc , for which we assume v esc = 550 km/s [38]. The solar circular speed is by convention taken as v 0 = 220 km/s and we boost from the galactic frame to the detector rest frame with v E = (0, v 0 , 0) + v pec + v e , where v pec = (11.1, 12.2, 7.3) km/s [39] and we use the expression for v e from [40].
The minimum speed to recoil with an energy E R additionally depends on the excitation energy E * : where µ A is the nucleus-dark matter reduced mass. The minimum speed is larger for bigger E * since part of the kinetic energy of the incoming dark matter particle is required to excite the nucleus. This means that for the same E R , elastic and inelastic scatter processes probe different parts of f DM ( v) [25].
In this paper we only consider axial-vector interactions of the type where χ is the dark matter (here assumed to be a fermion), ψ q are the light-quark fields (q = u, d, s) and A q are the (model-dependent) dark matter-quark coupling constants. The total spin-dependent differential cross-section applicable for this operator can generally be written as where µ n is the nucleon-dark matter reduced mass, the sum is over the isotopes that have spin, f A is the fractional abundance of the xenon isotope, J A is the ground-state spin of the isotope (J 129 = 1/2 and J 131 = 3/2) and σ 0 n is the elastic cross-section to scatter off a pointlike neutron in the limit of zero-momentum transfer. The structure factors S A (E R ) describe how the dark matter interacts with the nucleus and depend on the isotope. We take the central values of the one + two-body expressions from [24] and [25] for elastic and inelastic scattering respectively. In both cases we only consider the neutron structure factors S n A (E R ) since the proton structure factors are always at least a factor of 10 smaller.
The recoil spectra as a function of the xenon nucleus's recoil energy E R are shown in figure 2. The left and right panels show the spectra for m DM = 100 GeV and σ 0 n = 10 −40 cm 2 , and m DM = 1000 GeV and σ 0 n = 10 −39 cm 2 respectively. The total elastic and inelastic spectrum are shown by the black dotted and black dashed lines respectively. The orange �� /�� � [������/�/��/���] Figure 2. The left and right panels show the recoil spectra for the elastic and inelastic processes in terms of E R , the energy of the recoiling nucleus, for two values of the dark matter mass m DM and the scattering cross-section σ 0 n . The various curves show the individual rates for the xenon isotopes that participate in the scattering for the axial-vector (spin-dependent) interaction, namely 129 Xe and 131 Xe, together with the total rate. The elastic rate always dominates implying that the initial discovery of dark matter will always be made with the elastic scattering process. The difference between the elastic and inelastic rates are smaller for heavier particles suggesting that it should be easier to find evidence for the inelastic process with heavier particles. We remind the reader that E R is not directly measured in a two-phase xenon detector, rather, it is the scintillation signals S1 and S2. and green lines show the contribution of 129 Xe and 131 Xe to the elastic spectrum, while the blue and red lines show the contribution of 129 Xe and 131 Xe to the inelastic spectrum. The elastic spectrum is always larger than the inelastic spectrum with the most noticeable difference at small E R . The inelastic spectrum drops to zero at small E R because energy and momentum conservation do not allow for the xenon nucleus to be excited while remaining at rest after the dark matter interaction. The larger elastic scattering rate implies that for the axial-vector interaction, a discovery of dark matter will always first be made with the elastic scattering process.
The elastic spectrum is dominated by scattering with 129 Xe at low energy, while scattering with 129 Xe dominates for all energies in the inelastic spectrum. Comparing the left and right panels, we see that the inelastic spectra are closer to the elastic spectra for m DM = 1000 GeV. At low energies and for the mass values shown, the elastic spectra display the characteristic scaling dR/dE R ∝ σ 0 n m −1 DM . This scaling does not continue at higher recoil energies because for smaller masses, it is only the particles in the tail of the Maxwell-Boltzmann distribution that have sufficient kinetic energy to induce higher recoil energies of the xenon nucleus, thus producing an additional suppression in the rate (this is manifested mathematically through a higher value of v min for smaller m DM ). The inelastic spectra also do not show the characteristic scaling at any energy for these masses. This is for the same reason as in the elastic case, namely, many more incoming dark matter particles have a larger kinetic energy for higher masses. This is especially noticeable for the 131 Xe spectra where there is a factor ∼ 6 difference between the 129 Xe and 131 Xe spectra at m DM = 100 GeV, while only a factor ∼ 3 at m DM = 1000 GeV. This suggests that it should be easier to find evidence for the inelastic process for heavier dark matter particles.

JCAP05(2016)033
This discussion so far has only considered the recoil energy of the nucleus and has not accounted for the energy deposited by the photon from the de-excitation process. Although the de-excitation will not change the total (integrated) scattering rate, it must be accounted for when modelling the signal that a two-phase xenon detector measures. The next subsections address how we model this.

Generating light and charge signals
Two-phase xenon detectors do not directly measure the energy. Instead, a particle interacting in the fiducial volume of the liquid xenon produces two measurable signals referred to as the S1 and S2 signal. 1 An interaction in the liquid xenon produces ions and excitons which produce photons and electrons. The quantity S1 is a measure of the number of photoelectrons (PE) from the prompt scintillation due to the photons in the liquid xenon. The electrons are drifted in an electric field to the xenon gas phase, where they are extracted and accelerated. These extracted electrons create a secondary scintillation, denoted as S2. Electronic and nuclear events produce different characteristic S1 and S2 signals, which allows two-phase xenon detectors to discriminate between these two event classes. Canonical dark matter interactions result in nuclear events while most background events are electronic events. This ability to discriminate between nuclear and electronic events is one important reason why two-phase xenon detectors have been so successful in constraining the dark matter scattering cross-section.
For an energy deposition E, the expectation values for S1 and S2 can be expressed as S1 = g 1 n γ (E) and S2 = g 2 n e (E) respectively, where the measurement gains g 1 and g 2 relate the number of produced photons n γ (E) and electrons n e (E) to the expected number of detected PEs. 2 The gain g 1 is the probability that a photon produced at the centre of the detector strikes a photomultiplier tube (PMT) and is converted to a PE. The gain g 2 = × Y is the product of the probability of extracting an electron from the liquid to the gas ( ) and the amplification factor (Y ) converting a single ionisation electron to photoelectrons. The S2 signal measured from the bottom PMTs (S2 b ) is usually used because the light collection efficiency is more homogeneous on these PMTs [41]. We therefore use S2 b and assume that it is related to the total S2 signal by S2 b = 0.43×S2, as found in XENON100 and LUX [44,45]. We will discuss realistic values of g 1 and g 2 in the next subsection and for now, leave them as free parameters in our discussion.
To simulate signal processes observed by a two-phase xenon detector in a realistic fashion, we generate the signal with a Monte Carlo simulation along the lines of ref. [46]. We use the NEST phenomenological model [33][34][35][36] to model the average number of photons n γ (E) and electrons n e (E) produced by an electronic-or nuclear-type interaction. In addition to their dependence on the energy E, n γ (E) and n e (E) also depend on the electric drift field applied across the liquid, which varies for different experiments. The specific formulae used in our modelling are given in appendix A. We must include fluctuation effects, which can be divided into two types: intrinsic and detector fluctuations. We discuss the implementation of each in turn beginning with the intrinsic fluctuations.
Our signal generation begins by drawing the energy E of the incident particle from the input energy spectrum. For dark matter events, this is simply the recoil spectrum eq. (2.1)

JCAP05(2016)033
(as in figure 2), while the background distributions are discussed in section 3 (displayed in figure 5). For this energy, we find the total number of quanta N quanta by drawing from a Normal distribution with mean n quanta and variance F · n quanta , where F = 0.05 is the Fano factor [47]. We next separate N quanta into excitons and ions. The number of ions N i is drawn from a binomial distribution with N quanta trials and a probability (1 + n ex /n i ) −1 that an ion is produced. The number of excitons N ex is simply N ex = N quanta − N i . Our expressions for n γ (E) and n e (E) also include the effect of recombination fluctuations; we assume that the number of ions that recombine N recom i follows a Normal distribution with mean rN i and variance σ 2 R = (1 − r)CN 2 i , where C = 0.0056 [36]. Our final result is that n e (E) = N i − N recom i and n γ (E) = f l (N ex + N recom i ), where f l is a quenching factor. The Monte Carlo process is the same for both nuclear and electronic recoils; the difference is that n quanta , n ex , n i , r and f l differ for nuclear and electronic recoils. The calculation of the mean quantities n quanta , n ex , n i , r and the quenching factor f l are described in appendix A. There is also a small difference between gamma-and beta-electronic interactions that we account for by rescaling n γ (E) and n e (E) calculated for a beta-interaction to obtain the result for the gamma-interaction. We perform this rescaling at this point, after the intrinsic fluctuations. Further details are given in appendix A.
We next include detector fluctuations in our calculation of S1 and S2 b . For S1, the number of photoelectrons N PE is drawn from a binomial distribution with n γ (E) trials and success probability g 1 . The final result for S1 also accounts for the PMT resolution: S1 is drawn from a Normal distribution with mean N PE and variance σ 2 PMT N PE . For S2 b , the number of electrons N e that are extracted from the liquid to the gas phase follows a binomial distribution with n e (E) trials and success probability . To account for the amplification factor from converting ionisation electrons to photoelectrons, we draw from a Normal distribution with mean 0.43 · Y · N e and variance σ 2 PE b N e (the factor 0.43 is the factor that relates S2 and S2 b ).
When generating the S1 signal from the inelastic scattering process, we combine the number of photons from the nuclear recoil with the photons from the de-excitation gammaray with energy E * after the intrinsic fluctuations (for which the two processes are treated independently) but before including the detector fluctuations. An analogous procedure is performed for S2 b except it is electrons that we combine before including the detector fluctuations.

Two-phase xenon detector parameters
The LUX and XENON collaborations have produced the most sensitive xenon detectors. The current LUX detector has a fiducial mass of 118 kg and they have collected an exposure of 0.028 tonne-years [45], with an ultimate aim of collecting around 0.2 tonne-years [12]. The applied drift field of 181 V/cm is lower than in previous experiments. For instance, ZEPLIN-III [48] and XENON100 [49] had fields of 3400 V/cm and 530 V/cm respectively. However, the LUX light collection efficiency is much higher than in previous detectors, corresponding to a value of g 1 ≈ 0.12 PE/γ [50,51]. Unfortunately the electron extraction efficiency is lower than was anticipated, with ≈ 50% [50,51]. The amplification factor is Y ≈ 24.6 PE/e [45] so that g 2 ≈ 12 PE/e. The energy resolutions are σ PMT ≈ 0.5 PE/γ and σ PE b ≈ 3.6 PE/e [52].
XENON1T is the successor to XENON100. It will have a fiducial mass of approximately 1 tonne, a design drift field of 1000 V/cm and a light collection efficiency similar to LUX [13]. It is expected that the extraction efficiency will be 100%, as achieved in XENON100. The amplification factor and resolutions will be similar to those in LUX and XENON100 [53].

JCAP05(2016)033
The follow-up to LUX is LZ [14,15], with a projected fiducial mass of approximately 5.6 tonnes and a drift field of 700 V/cm [54]. XENONnT is the successor of XENON1T and is designed to have similar characteristics as XENON1T but with a total mass of approximately 7 tonnes [16]. As XENONnT and LZ will run for a number of years, an exposure of 15 tonneyears is readily achievable. Finally, there are plans for DARWIN, a much larger experiment with a fiducial mass of around 20 tonnes [18]. Studies assuming a drift field of 500 V/cm and an ultimate exposure of 200 tonne-years have been performed [43]. This large exposure gives an indication of the ultimate reach of xenon detectors.
Future collaborations will obviously aim to optimise their respective detectors. It may be difficult to increase or even maintain the light collection efficiency because larger detectors collect a smaller fraction of the scintillation signal. LZ's proposal is that g 1 > 0.075 [15] and DARWIN studies have assumed the value reached in LUX [43]. It should be possible to maintain an extraction efficiency close to unity and the amplification factor may be as large as Y = 50 PE/e [15].
In the remainder of this paper, we show results for two benchmark scenarios, XenonA200 and XenonB1000, which should bracket the expected performance of upcoming experiments: • XenonA200 corresponds to a detector with lower drift field and lower extraction efficiency. We assume a drift field of 200 V/cm and the parameters g 1 = 0.07 PE/γ, = 50%, Y = 25 PE/e so that g 2 = 12.5 PE/e.
• XenonB1000 corresponds to a detector with a higher drift field and perfect extraction efficiency. We assume a drift field of 1000 V/cm and the parameters g 1 = 0.12 PE/γ, = 100%, Y = 50 PE/e so that g 2 = 50 PE/e.
In both scenarios, we assume that σ PMT = 0.5 PE/γ and σ PE b = 3.6 PE/e and that S2 b = 0.43 × S2. Our benchmark exposure is 15 tonne-years unless stated otherwise. Finally, we assume that all measurement efficiencies are 100% since the signals of interest are far from thresholds (where the efficiencies begin to deviate from 100%).

Observable signals and their rate
Having described our procedure for generating light and charge signals, we can proceed to generate the observable signals for our two benchmark scenarios. The solid and dashed contours in figure 3 show where 68% and 95% of events occur in the S1 -S2 b plane for fixed input energies. The left and right panels correspond to the XenonA200 and XenonB1000 benchmark scenarios, respectively. The keV β ER , keV γ ER and keV NR labels in figure 3 indicate that the energy originates from a beta-electronic, gamma-electronic and nuclear recoil, respectively. The black contours show the signal region for a nuclear recoil of 40 keV, the red and blue contours show the signal region for a 39.6 keV electronic event induced by a beta-and gamma-ray respectively, and the purple and pink contours show the signal region for a 80.2 keV electronic event induced by a beta-and gamma-ray, respectively. The brown and orange contours show the signal region for an event that occurs for inelastic scattering: in this case the nuclear recoil energy is 40 keV and the gamma-electronic energy is 39.6 keV and 80.2 keV, corresponding to the energies of the photon emitted in the de-excitation of the 129 Xe and 131 Xe isotopes, respectively.
We first discuss the features common to both panels. These features are well known properties of two-phase xenon detectors and are reviewed in much more detail in [55]. We see that a nuclear recoil typically produces a much smaller S1 and S2 b signal compared Figure 3. The solid and dashed contours show where 68% and 95% of events occur in terms of the observable scintillation signals S1 and S2 b for the fixed input energies indicated. The left and right panels show results for the XenonA200 and XenonB1000 benchmark scenarios. The keV β ER , keV γ ER and keV NR labels indicate that the energy originates from a beta-electronic, gamma-electronic and nuclear recoil. The brown and orange contours show the signal region for an event that occurs for inelastic scattering, which includes energy from both a nuclear recoil and a de-excitation photon. Inelastic signal events have higher S1 and S2 values than for elastic scattering events, for which the search window is typically S1 ≤ 30 PE. All of the contours are tilted because of recombination fluctuations. Note the change of scale on both axes between the two panels.
to an electronic recoil of the same energy. The usual XENON100 and LUX dark matter searches for elastic scattering define a S1 search window up to 30 PE [45,49]; we see that for inelastic signals, we will have to consider much higher values of S1. The difference between a gamma-and beta-interaction of the same energy is relatively small, O(10%), for both S1 and S2 b . It is also apparent that adding a nuclear recoil to an electronic recoil only slightly increases the S1 and S2 b signals compared to a pure electronic recoil. Both panels show that the contours are tilted, matching the behaviour observed with real data (see e.g. [56]). This is especially obvious in the events where the S1 and S2 b signals are dominated by electronic recoils while for nuclear recoils, the tilt is much smaller. The origin of the tilt is well-known: it is a result of recombination fluctuations, which are 100% anti-correlated in scintillation (S1) and ionisation (S2) [57]. In contrast, the detector fluctuations smear along constant S1 and constant S2 only. The tilt is smaller for the nuclear recoil region because the detector fluctuations are larger than recombination fluctuations.
We next discuss the features that differ between the panels. The first important difference is that the S1 and S2 b signals are much larger in the right panel for all configurations (note that the scales differs in the two panels). The larger S2 b signal is a result of two effects. The first is that the extraction efficiency and amplification factor Y are both twice and therefore g 2 is four times larger for XenonB1000. If this were the only effect, then S2 b would be four times larger for the same input energy. In fact, we see that S2 b is around six times larger in the right panel. This is because the larger drift field also increases S2 b . The larger drift field reduces the recombination fraction so more of the ions survive to form the S2 b signal. The higher drift field also reduces the S1 signal for the same reason; fewer ions recombine producing fewer prompt scintillation light. However, g 1 is 70% larger for XenonB1000, which compensates for the reduction from the higher drift field. This is why Figure 4. The left and right panels show the recoil spectra for the elastic and inelastic processes in terms of S1 for the two benchmark scenarios XenonA200 and XenonB1000 described in section 2.3. These spectra are for m DM = 100 GeV and σ 0 n = 10 −40 cm 2 . While the elastic spectrum falls off rapidly, the inelastic spectrum has two distinct peaks. The first peak is dominated by the 39.6 keV de-excitation from 129 Xe while the second peak is dominated by the 80.2 keV de-excitation from 131 Xe. Similar to the left panel of figure 2, the inelastic rate is suppressed by about two (three) orders of magnitude with respect to the elastic rate for the 129 Xe 131 Xe isotope.
the S1 values are actually ∼ 20% larger in the right panel. Finally, an important difference is that there is less overlap between the contours from an inelastic signal (brown and orange contours) and the contours from a potential beta-background source (red and purple contours) for XenonB1000. This greater separation is again an effect of the higher drift field. We will see in section 5 that this better discrimination is ultimately responsible for the greater sensitivity of the XenonB1000 benchmark scenario.
We previously only gave the differential event rate in terms of the xenon recoil energy E R (cf. eq. (2.1) and figure 2). We are now in a position to calculate the event rate in terms of the observable quantities S1 and S2 b . The differential event rate is [41] where dR/dE R is eq. (2.1) and pdf(S1, S2 b |E R ) is the probability density function, which we generate with the Monte Carlo process described in section 2.2.
As an example of our results, we show in figure 4 the differential event rate in terms of S1 for m DM = 100 GeV and σ 0 n = 10 −40 cm 2 (additionally, the results for m DM = 1000 GeV and σ 0 n = 10 −39 cm 2 are shown in figure 6). The left and right panels correspond to the XenonA200 and XenonB1000 benchmark scenarios, respectively. The dR/dS1 spectrum is obtained by additionally integrating eq. (2.5) over S2 b . The colour scheme of the lines matches figure 2: the total elastic and inelastic spectrum are shown by the black dotted and black dashed lines respectively. The orange and green lines show the contribution of 129 Xe and 131 Xe to the elastic spectrum, while the blue and red lines show the contribution of 129 Xe and 131 Xe to the inelastic spectrum. As in the left panel of figure 2, the inelastic rate for scattering with 129 Xe 131 Xe is suppressed by about two (three) orders of magnitude with respect to the elastic rate.  Figure 4 shows the well known fact that the elastic spectrum falls off rapidly with S1. In contrast, the inelastic spectrum has two distinct peaks whose origin is clear. The first peak is due to the 39.6 keV de-excitation photon from the 129 Xe isotope while the second peak is from the 80.2 keV de-excitation photon from the 131 Xe isotope. The peak S1 values agree with the values shown in figure 3. The peak differential rate is slightly higher for the inelastic process in the left panel (corresponding to XenonA200 ) because the spectrum is slightly more peaked in S1 (the integrated rate is the same).

Background rates
Our ultimate aim is to assess the discovery potential of the inelastic signal. In order to do this, the background signals must be quantified. A comprehensive study of the backgrounds for tonne-scale xenon detectors was performed in [58] and similar rates and distributions were also presented by the LZ collaboration [15,59]. We summarise the relevant results for our study and refer the reader to [15,58,59] for further details.
The dominant backgrounds are those that produce electronic recoils in the signal range of interest, S1 ≤ 600 PE, which corresponds to an energy range of approximately 0-300 keV. As discussed in [58,59], the background rates that dominate in order of decreasing importance are the 2νββ-decay of 136 Xe, elastic neutrino-electron scattering from pp and 7 Be solar neutrinos, decays of 85 Kr and 222 Rn and finally, radioactivity from detector materials. All of these backgrounds are beta-electronic sources [23]. The background rates and their uncertainties used in this study are shown in figure 5. We comment on each of the rates in turn.
We recalculated the background rates from the 2νββ-decay of 136 Xe and the elastic scattering from pp and 7 Be solar neutrinos using updated parameters. For the 2νββ-decay, the neutrinos escape the detector while the beta-particles contribute to the background rate. We calculated the rate assuming the 136 Xe abundance is that of natural xenon: 8.86%. We use the most accurate measurement of the 136 Xe half-life by EXO-200, who found -11 -JCAP05(2016)033 [60], where we only quote the dominant systematic error. We use the distribution of the summed energies of the beta-particles from [61]. Our 2νββ rate is in good agreement with that shown in [58,59]. The dominant uncertainty of this rate is from T 1/2 , at the level of 3%. In comparison, the abundance of the 136 Xe isotope can be measured with a 0.05% accuracy [62].
The largest flux of solar neutrinos is from the pp chain. The pp flux measured by Borexino, (6.6 ± 0.7) × 10 10 cm −2 s −1 [63], is in good agreement with the prediction of the Standard Solar Model (SSM) 6.03 × (1 ± 0.06) × 10 10 cm −2 s −1 [64]. The SSM prediction is well understood so we use the theoretical flux and error in our calculation. The second largest rate from solar neutrinos is from 7 Be neutrinos. Borexino measured a flux of (4.84 ± 0.24) × 10 9 cm −2 s −1 [65], also in good agreement with the SSM prediction [64]. We use the measured value and error in our calculation. We use the neutrino-electron scattering crosssection from [66]. Finally, we use the electron neutrino survival probabilities listed in [67]. Our pp spectrum agrees well with [58,59]. Our 7 Be spectrum agrees well with [59] but is about a factor 1.6 smaller than in [58]. The origin of the discrepancy is unclear. 3 The third largest rate from solar neutrinos is from 13 N neutrinos. This rate is approximately 300 times smaller than the pp rate so it is a good approximation to ignore the contribution to the rate from all solar neutrinos except the pp and 7 Be neutrinos.
We use the 85 Kr, 222 Rn and detector material background rates from the study in [58], which assumes a 85 Kr contamination of 0.1 ppt and a 222 Rn level of 0.1 µBq/kg. While the 222 Rn rate is similar in the LZ study, the 85 Kr rate is a factor four smaller [68]. We use the result from [58] because the assumptions entering the calculation are clearer. XENON100 and EXO-200 have measured the 85 Kr contamination and 222 Rn level with an accuracy of 17% [69] and 10% [60] respectively, and we assume the same accuracy will be achieved in the future.
The detector material background rate is reduced by self-shielding of the liquid xenon so larger detectors, which have more xenon with which to shield, have a smaller rate. The rate reported here was for a DARWIN study and assumes a 14 tonne fiducial mass. The rate for LZ with a 5.6 tonne fiducial mass is about three times larger [59]. As before, we use the result from [58] because the assumptions entering the calculation are clearer. In any case, as this rate is always subdominant, a factor three difference has an almost negligible impact on our results. Both the material and 222 Rn background rates begin to increase after ∼ 170 keV however they always remain subdominant to the 2νββ-decay rate [58]. The detector material rate for XENON1T is predicted to 10% [68] and we assume the same accuracy will be achieved in future experiments.
Before leaving this sub-section, we briefly comment on background sources that do not produce electronic recoils. It is also possible that neutrons may directly excite the 129 Xe and 131 Xe isotopes creating another background source (that neutrons can excite the signal is actually an advantage since at least in principle, it allows the signal region to be calibrated in a real detector). The self-shielding of the liquid xenon and a dedicated muonveto system to reject muon-induced neutrons means that the neutron rate can be reduced to less than 5 × 10 −5 counts/t/yr/keV for single-scatter neutrons that elastically scatter of xenon [43]. For comparison, the dark matter signal rate that we consider in this paper is ∼ 10 −2 counts/t/yr/keV for σ 0 n 10 −40 cm 2 (cf. figure 6). Although there are no detailed studies that discuss inelastic neutron scattering (and such a study is beyond the scope of this paper), the inelastic neutron scattering cross-section is generally of the same order of  The dominant background rate is from the 2νββ-decay of the 136 Xe isotope (solid purple line), which is always at least 30 times larger than the dark matter signal. These panels demonstrate that observing the inelastic signal with a single-phase xenon experiment that only measures the S1 signal will be very challenging because of the large background rate. magnitude as the elastic scattering cross-section [70]. Therefore, we assume that the inelastic neutron scattering rate is similar to the elastic scattering rate and is therefore always significantly smaller than the dark matter signal rate, so we ignore this background contribution in our study. A more detailed study to confirm this assumption is desirable.

Comparing background and signal rates
We now have everything to model the signal and the background for the XenonA200 and XenonB1000 detector scenarios. The left and right panels in figure 6 show the background and signal rates as a function of S1 for the two benchmark scenarios. The red and blue lines show the signal rate corresponding to inelastic scattering off the 129 Xe and 131 Xe isotopes for m DM = 1000 GeV and σ 0 n = 10 −39 cm 2 . Comparing the background rates in the left and right panels, we see only minor differences. In both cases the dominant background is from the 2νββ-decay of 136 Xe (solid purple line). The most obvious difference is in the rate from detector materials (dotted grey line), where for XenonB1000, the rate is higher for S1 values corresponding to the 80.2 keV de-excitation. The main point to take away from both panels of figure 6 is that the signal rate is always at least 30 times smaller than the background rate. This demonstrates that observing this signal with a single-phase xenon experiment that only measures the S1 signal will be very challenging.
Two-phase experiments provide additional information in the form of the S2 signal. In figure 7 we therefore plot the signal and background distributions in the log 10 (S2 b /S1)-S1 plane traditionally used by two-phase xenon experiments. The black and purple lines show the electronic and nuclear recoil bands, respectively. The solid lines show the median while the dashed lines show ±1.28σ around the median, such that 10% of events are above and 10% below the dashed lines. The bands are calculated by passing a constant energy spectrum through our detector simulations for nuclear and beta-electronic recoils. The overall shape of -  the bands, and in particular that they separate at large S1, matches the behaviour observed with real detectors (see e.g. [46]). The blue and red contours indicate where 68% and 95% of events occur for inelastic scattering off the 129 Xe and 131 Xe isotopes for m DM = 1000 GeV and σ 0 n = 10 −39 cm 2 , respectively. Unlike the contour regions shown in figure 3, these contours are not elliptical but have a more extended shape. This shape change arises because figure 7 includes the effect of all possible recoil energies of the nucleus while figure 3 was for a single nuclear recoil energy. Finally, the circles and triangles show the simulated events expected for an exposure of one tonne-year and the dark matter parameters mentioned above. The open grey circles show the electronic background events, which as expected from figure 6, become more abundant at higher values of S1. The filled blue and filled red triangles show the inelastic events from the 39.6 keV and 80.2 keV de-excitation after scattering off the 129 Xe and 131 Xe isotopes. The open green triangles show the events from elastic scattering off xenon, which are more abundant at smaller values of S1.
Both panels of figure 7 show that the signal and background distributions are slightly displaced. The displacement occurs for two reasons. The first is that nuclear recoils also have a lower S2 b for the same S1 compared to an electronic event (this is why the nuclear band is below the electronic recoil band) and the second is that gamma-electronic interactions have a higher S1 and lower S2 b than a beta-interaction of the same energy (as shown in figure 3). Both effects mean that the inelastic signal region lies below the electronic recoil band. This displacement is crucial as it allows for some discrimination between signal and background events. This means that two-phase xenon detectors should have a significantly better sensitivity than single-phase detectors.

JCAP05(2016)033
Finally, we discuss the differences between the two benchmark scenarios. All of the signals have larger S1 and S2 b values for the XenonB1000 scenario because of the larger g 1 and g 2 values. The extent to which the signal regions extend below the electronic recoil band depends on the detector parameters, particularly the applied electric drift field. For XenonA200 (left panel) which has a drift field of 200 V/cm, 78% of the 129 Xe inelastic events signal fall below the lower dashed line of the electronic recoil band. In comparison, for XenonB1000 (right panel) where the drift field is 1000 V/cm, 92% of the 129 Xe inelastic events signal fall below this line. For the 131 Xe signal region, 92% of the events fall below the lower dashed line of the electronic recoil band for both the XenonA200 and XenonB1000 benchmark scenarios so we expect a similar sensitivity to the inelastic signal from the 131 Xe isotope for these scenarios.
To demonstrate that the drift field is primarily responsible for the separation of the signal region and electronic recoil band, we repeated this analysis with the detector parameters of XenonA200 but with a drift field of 1000 V/cm instead of 200 V/cm. In this case, we found that 91% of the 129 Xe inelastic events signal fall below the lower dashed line of the electronic recoil band, similar to the 92% obtained for XenonB1000. Similarly, for the detector parameters of XenonB1000 but with a drift field of 200 V/cm instead of 1000 V/cm, we found a value of 80%, similar to the value 78% obtained for XenonA200. Figure 7 was generated for m DM = 1000 GeV and σ 0 n = 10 −39 cm 2 but similar signal regions hold for other masses. This should not be too surprising since most of the S1 and S2 b signal originates from the de-excitation photon whose energy is always the same. The primary change is that the ratio of the 129 Xe to 131 Xe rate is larger for smaller mass values (cf. figures 4 and 6).

Characterising the detection sensitivity
In this section we describe our method for characterising the sensitivity of two-phase xenon experiments to the inelastic scattering process. We will do this by calculating the 'discovery limit' or as we will call it, the discovery reach. This was introduced in [71] and has been used extensively to characterise the limiting effect of the neutrino background (see e.g. [19]). We first describe the formalism behind this frequentist approach and then provide specific details of our calculation.
The discovery reach is the smallest cross-section for which 90% of experiments make at least a 3σ discovery of the signal under consideration. To calculate it, we make use of the frequentist test statistic for the discovery of a positive signal [72]: where the profile likelihood ratio is and the hats (ˆ,ˆ) indicate that the parameters are those that maximise the extended likelihood L. Here A BG = {A 2νββ , A pp , A Kr , A Rn , A Be , A mat } are the amplitudes of the six background components discussed in section 3.

JCAP05(2016)033
In our case the extended likelihood [73] (for a given value of the dark matter mass) is where µ DM ∝ σ 0 n and µ BGj ∝ A BGj are the mean number of events from dark matter and the background processes respectively, f DM and f BG are the unit normalised twodimensional probability distribution functions for the signal and background processes in the S1-log 10 (S2 b /S1) plane, N is the total number of observed events and {S1 i , log 10 (S2 b /S1) i } are the values for a single event. Finally, L m (A BGm ) are the individual likelihood functions for the background normalisations, which we assume are Normal distributions with a standard deviation given by the respective error quoted in figure 5. As we are dealing with hypothetical experiments, we generate the unit normalised two-dimensional probability distribution functions f DM and f BGj from Monte Carlo by generating approximately two million events for each process in the S1-log 10 (S2 b /S1) plane.
The results of Wilks [74] and Wald [75] allow us to relate the significance with which we can reject the background-only hypothesis (σ 0 n = 0) to the test statistic in a simple way: where Z 0 is the number of standard deviations. In appendix B, we explicitly demonstrate that the approximation of Wald is good so that eq. (4.4) is accurate. To obtain the discovery reach for each value of m DM , we simulate a minimum of 2500 mock experiments and find the cross-section σ 0 n for which 90% of experiments have Z 0 ≥ 3. We will present a separate discovery reach for inelastic scattering off the 129 Xe and 131 Xe isotopes. We are able to do this because, as figure 7 shows, the two signal regions are well separated from each other and also from the nuclear recoil band, which contains events from the elastic scattering process. The profile likelihood analysis takes into account the expected dark matter signal in the S1-log 10 (S2 b /S1) plane so no cuts to identify a signal region are required. However in practice, to improve the run-time of our calculations, for each discovery reach calculation we restrict our analysis to the S1 and log 10 (S2 b /S1) values around the dark matter signal region of interest. The restricted region is chosen to contain at least 95% of the events for each dark matter signal under consideration. By trying different regions, we found that our results are not sensitive to any reasonable choice. In appendix C, we also provide a discovery reach calculation using a more conservative cut-and-count method. This serves as a useful cross-check against the profile likelihood analysis.

Discovery reach for two-phase xenon detectors
We present in figure 8 the main results of this paper. The red and blue lines show the discovery reach for detecting inelastic scattering off the 129 Xe and 131 Xe isotopes for an exposure of 15 tonne-years. These lines show the smallest cross-section for which 90% of experiments are able to make at least a 3σ discovery of the inelastic signal. The left and right panels Figure 8. This figure shows the sensitivity of two-phase xenon experiments to the inelastic scattering process, which is our main result. The blue and red lines in both panels show the discovery reach for inelastically scattering off the 129 Xe and 131 Xe isotopes, respectively. The left and right panels show the results for the XenonA200 and XenonB1000 benchmark scenarios assuming a 15 tonne-year exposure. In the parameter space above these lines, 90% of experiments will make at least a 3σ detection of the inelastic signal. The discovery reach for inelastically scattering off the 129 Xe isotope is better for the XenonB1000 scenario, while for scattering off 131 Xe, both benchmark scenarios have similar sensitivity. Also shown is the LUX exclusion limit (black dashed line) from their search for elastically scattering dark matter and the projected exclusion limit from XENON1T assuming a two tonne-year exposure (black dot-dashed line).
show the results for the XenonA200 and XenonB1000 benchmark scenarios (described in section 2.3). The black dashed line shows the LUX 90% CL limit on the spin-dependent dark matter-neutron cross-section from their reanalysis of the 2013 search for dark matter that elastically scatters with the xenon isotopes [45,76,77]. The black dot-dashed line shows the projected limit from the XENON1T search for elastically scattering dark matter particles after an exposure of two tonne-years, which should be achieved by 2018.
For both xenon isotopes and both benchmark scenarios, figure 8 shows that the discovery reach of the inelastic signal is below the current LUX exclusion limit for a dark matter mass greater than ∼ 100 GeV. This means that for dark matter particles that are heavier than this, it is possible for the inelastic signal to be detected by a future two-phase xenon detector that collects an exposure of 15 tonne-year (such as LZ or XENONnT). The parameter space where the inelastic signal may be detected is populated by many dark matter models, including neutralino scenarios where the higgsino component is large (see e.g. [78][79][80]). XENON1T is expected to be significantly more sensitive than LUX and will probe all of the parameter space where the inelastic signal may be detected with a 15 tonne-year exposure. Therefore, if the inelastic signal is ever to be detected with this exposure, XENON1T should find evidence for the elastic scattering dark matter signal by 2018.
Comparing the left and right panels of figure 8, we see that the discovery reach for inelastic scattering off the 131 Xe isotope (red line) is similar for both benchmark scenarios. This is because the ability to discriminate between signal and background processes is similar for both scenarios. In contrast, the discovery reach for inelastic scattering off the 129 Xe isotope (blue line) is a factor ∼ 3.5 lower for the XenonB1000 benchmark scenario. This is because, as figure 7 shows, more of the signal region lies below the electronic recoil band for  the XenonB1000 scenario so the discrimination power is better (cf. the discussion at the end of section 3.1 where the discrimination power was quantified). The discovery reach shown in figure 8 assumes an exposure of 15 tonne-years and the backgrounds rates discussed in section 3. We now explore how the discovery reach changes as we vary these assumptions. Firstly, we examine variations in the background rate. As we showed in section 3, the dominant background is the from the 2νββ-decay of the 136 Xe isotope so it is possible to reduce this background rate by reducing the abundance of the 136 Xe isotope. Depleting (or enriching) the 136 Xe isotope from xenon is relatively straightforward as demonstrated by experiments that use xenon enriched in 136 Xe to search for neutrinoless double beta decay. The upper two panels of figure 9 show how the discovery reach changes as we vary the 136 Xe abundance. The left and right panels show the results for the XenonA200 and XenonB1000 benchmark scenarios, respectively. The results are shown for two dark matter mass values and we plot the discovery reach cross-section normalised to the discovery reach assuming that the fractional abundance of 136 Xe is 8.86%, which is the abundance in -18 -

JCAP05(2016)033
natural xenon and is the value that we have assumed throughout the paper. As expected, the discovery reach extends to smaller values of the cross-section as the fractional abundance is reduced. Figure 9 shows that the variation does not depend on the dark matter mass and is only weakly dependent on the benchmark scenario. For both scenarios, we see that lowering the 136 Xe abundance to 1% means that the smallest cross-section for which the inelastic signal may be discovered is reduced by ∼ 35%.
Secondly, we examine variations in the exposure. These results are shown in the lower two panels of figure 9, where the discovery reach cross-section has been normalised to the discovery reach assuming a 15 tonne-year exposure. Figure 9 again demonstrates that the variation does not depend on the dark matter mass and is only weakly dependent on the benchmark scenario. The discovery reach for inelastically scattering off 129 Xe for the XenonA200 scenario decreases more slowly as the exposure increases compared to the other scenarios because this signal region is most dominated by background processes (cf. the discussion at the end of section 3.1). For a background free signal region, the improvement in the discovery reach is expected to scale as (exposure) −1 . As there is always some background contamination for the inelastic signals and the benchmark scenarios that we consider, we instead find that the discovery reach scales as ∼ (exposure) −0.7 . In practice, this means that for a 200 tonne-year exposure, a benchmark exposure used in sensitivity studies for DARWIN, the various discovery reach lines presented in figure 8 should be lowered by a factor ∼ 5.
We end this section by returning to the question of whether a two tonne-year exposure of XENON1T can probe all of the parameter space where the inelastic signal may be discovered. With the scaling of the exposure determined in figure 9, we find that the exposure required to reach the XENON1T exclusion limit for all scenarios except the 129 Xe signal in the XenonB1000 benchmark scenario is ∼ 500 tonne-year. Such a large exposure is unlikely to be achieved in the foreseeable future. For the 129 Xe signal in the XenonB1000 scenario, an exposure of approximately {225, 90, 70} tonne-year for m DM = {150, 1000, 10000} GeV is required for the discovery cross-section to reach the XENON1T limit shown in figure 8. The exposure can be reduced to approximately {165, 60, 45} tonne-year if the 136 Xe abundance is reduced to 1%. This demonstrates that it may be possible to discover inelastic scattering off the 129 Xe isotope even for cross-sections below the XENON1T limit, but only for optimal detector parameters (as in the XenonB1000 scenario) and with large exposures that will only be achieved with detectors such as DARWIN.

Conclusions and outlook
The canonical search for dark matter with direct detection experiments is for an elastic scattering process where the dark matter simply causes the nucleus to recoil. It was long ago realised that low-lying inelastic transitions of the nucleus may also play an important role since the dark matter's kinetic energy is sufficient to excite the target nucleus. In this instance, rather than just measuring the recoil of the nucleus, direct detection experiments measure the nuclear recoil energy together with the photon-energy released when the nucleus transitions back to the ground state (see figure 1).
The inelastic scattering rate does not have the nucleon-number-squared enhancement (∼ 10 4 ) found with elastic spin-independent interactions so the inelastic signal will only be measurable for spin-dependent interactions, whose elastic scattering rate also does not have the nucleon-number-squared enhancement. Two-phase xenon detectors are an excellent probe of the elastic and inelastic spin-dependent interaction having two isotopes sensitive to these -19 -

JCAP05(2016)033
processes, 129 Xe and 131 Xe, that each comprise approximately 25% of natural xenon and have low-lying excitations at 39.6 keV and 80.2 keV, respectively. The purpose of this paper was to quantify the sensitivity of future tonne-scale two-phase xenon experiments, such as LZ, XENONnT and DARWIN, to the inelastic signal. We do this for the axial-vector interaction (eq. (2.3)), for which accurate calculations of the nuclear structure functions are available.
We considered two benchmark scenarios, XenonA200 and XenonB1000 (described in section 2.3), whose most important difference is the applied drift field of 200 V/cm and 1000 V/cm, respectively. The parameters in these scenarios were chosen because they should bracket the performance of future experiments. We implemented a realistic Monte Carlo simulation of a two-phase xenon detector to model these scenarios, relying on the NEST phenomenological model to describe the interactions of the nucleus and photon in liquid xenon. This was vital so that we could translate energies into the measurable quantities, the primary (S1) and secondary (S2) scintillation signals (see figure 3). We also had to quantify the background rates, finding that the 2νββ-decay of 136 Xe dominates (see figure 5).
We demonstrated that two-phase xenon detectors allow for some discrimination between the inelastic signal and the background events because the signal region has a smaller log 10 (S2 b /S1) value compared to the main backgrounds (see figure 7). Our main results were shown in figures 8 and 9, where we quantified the sensitivity of our benchmark scenarios to the inelastic signal in terms of the discovery reach, which is the smallest cross-section for which 90% of experiments detect the signal with at least a 3σ significance. This cross-section is below the current LUX exclusion limit for a dark matter mass greater than ∼ 100 GeV, implying that for dark matter particles that are heavier than this, it is possible for the inelastic signal to be detected with a future two-phase xenon detector. Except in the case of optimal detector parameters (as in the XenonB1000 scenario) and large exposures (more than 50 tonne-years), XENON1T, with a two tonne-year exposure, will probe all of the parameter space where the inelastic signal may be detected with their search for elastically scattering dark matter.
We end by discussing some of the possible extensions of this work. Firstly, we were only able to consider the inelastic signal from the axial-vector interaction since this is the only spin-dependent interaction for which the inelastic structure functions have been calculated. It would be desirable to calculate the discovery reach of other spin-dependent operators such as VA (χγ µ χψ q γ µ γ 5 ψ q ) or SP (χχψ q γ 5 ψ q ). Secondly, in this work we assumed that nuclear recoil and photon scintillation signals could not be distinguished. It may be possible that for some fraction of the events, the photon travels sufficiently far from the initial interaction to give a distinctive pulse shape different from background events. This would further improve the sensitivity to inelastic scattering signals. Thirdly, we focussed solely on two-phase xenon experiments as these are better able to distinguish between signal and background signals compared to single-phase xenon detectors. However it may be possible that tonne-scale single-phase detectors can improve their sensitivity to the inelastic signal by performing an annual modulation search or by modelling the S1 pulse shape. Finally, we have stated that a detection of the inelastic signal together with the elastic signal would point strongly to a spindependent interaction over a spin-independent interaction from a single xenon experiment. It would be desirable to have a dedicated study to quantify this statement and to concretely demonstrate how a measurement of both signals would help pin down the nature of the interaction between dark matter and the Standard Model particles.
-20 -I am particularly grateful to Nassim Bozorgnia for discussions that reignited my interest in this topic and for providing data from the EAGLE simulation, and to Andrew Brown for answering many of my naive questions. I'm also grateful for discussions with Rafael Lang and Simon Fiorucci at the Rencontres de Blois conference, with Martin Hoferichter, Philipp Klos and Achim Schwenk at the Mainz Institute for Theoretical Physics (MITP), with Alastair Currie and Marc Schumann at the TAUP conference, to Felix Kahlhoefer for comments on an early draft of this manuscript, and to Sebastian Liem for discussions on statistics. This work is part of the research programme of the Foundation for Fundamental Research on Matter (FOM), which is part of the Netherlands Organisation for Scientific Research (NWO).

A Mean photon and electron yields
We use NEST's semi-empirical model to calculate the mean number of photons and electrons from nuclear and electromagnetic interactions. In order that our results can be easily reproduced, in this appendix we collect the parameters and formulae that we use. The energy is denoted E and has units keV while the applied electric field is denoted by F and has units V/cm.

A.1 Nuclear recoils
Our treatment of nuclear recoils follows [36]. For an energy deposit E, the number of quanta n quanta is where n quanta = n i + n ex is the total number of ions n i and excitons n ex respectively and W = 13.7 eV is the average energy to produce a quanta. The quenching factor (to account for energy lost to atomic motion) is where k = 0.1394, = 11.5EZ −7/3 (Z = 54 for xenon) and g( ) = 3 0.15 + 0.7 0.6 + . The number of ions and excitons follows the ratio where ζ = 0.0472, α = 1.240 and β = 239. The probability that an ion recombines is where = γF −δ , γ = 0.01385 and δ = 0.0620. The number of electrons n e and photons n γ is given by is another quenching factor (accounting for the Penning effects, when two excitons interact to produce one exciton and one photon), η = 3.3 and λ = 1.14.

A.2 Electronic recoils from incident beta particles
The equivalent formulae for electronic recoils are generally simpler since there are no quenching factors. Our treatment follows [81]. The number of quanta is where W is the same as in eq. (A.1), the ratio of excitons to ions is the probability that an ion recombines is .10) and the number of electrons n e and photons n γ is The expression for˜ is more complicated in this case, depending on both E and F . We have that˜

A.3 Electronic recoils from incident gamma particles
There is a small difference between the electron and photon yields from incident gamma and beta particles [33]. The difference arises because an incident gamma produces multiple lower energy electron recoils. For instance, a 39.6 keV gamma typically results in Auger and photo-electrons with energies 4.5 keV, 4.8 keV, 5.3 keV and 25 keV respectively [46]. A complete description of these events requires a full simulation with NEST, which tracks the energy deposition of the incident gamma. Such a simulation is beyond the scope of this work. Instead we use the mean yields from a NEST simulation presented in figure 1 of [34] to rescale the yields from an incident beta particle: Here n γ [β(E, F )] is the yield from a beta-particle, given by eq. (A.12). We have made it explicit that the yields depend on both the incident energy E and drift field F . This rescaling is sufficient to also calculate the change in the number of electrons. Quanta conservation tells us that n e = n quanta − n γ (this follows from the usual assumption that energy lost to heat can be ignored for electromagnetic interactions), so increasing n γ for incident gamma particles automatically decreases n e , as shown in [34]. We only require the rescaling α factors at the de-excitation energies: 39.6 keV and 80.2 keV. In table 1 we give the factors at these energies and at the two values of the drift field that we consider in this study.  Table 1. The rescaling α γ factor to convert between the yields for incident beta particles and incident gamma particles. We give the values for the two de-excitation energies and for the values of the drift fields that we consider in this study. Figure 10. A comparison of the distribution of the test statistic q 0 under the background-only hypothesis from a Monte Carlo simulation (blue histogram) and from the asymptotic formula that follows from Wald and Wilks (black dashed line). The four panels show the result for the four signal regions that we consider: the signal regions from scattering off the 129 Xe and 131 Xe isotopes for the XenonA200 and XenonB1000 benchmark scenarios. The Monte Carlo and asymptotic formula agree well, with the implication that our use of eq. (4.4) is robust.

B Testing the approximations of Wilks and Wald
As discussed in [72], if Wald's approximation is good, then following the result of Wilks, the distribution of the test statistic q 0 under the background-only hypothesis H back should asymptotically follow In figure 10 we demonstrate that Wald's approximation is good by comparing f (q 0 |H back ) from 100000 Monte Carlo simulations (blue histogram) with eq. (B.1) (black dashed line). We find good agreement between the asymptotic distribution and our Monte Carlo distribution for the four signal regions of interest. This implies that our use of eq. (4.4) is justified.

C A cut-and-count analysis
In the calculation of the discovery limits in the main body of the paper, we employed a profile likelihood analysis that takes into account the position of each event in the S1-log 10 (S2 b /S1) plane. In this appendix, we present an alternative calculation of the discovery limit using a more conservative cut-and-count analysis. In this case, we define a signal box in the S1log 10 (S2 b /S1) plane and simply count the number of events that fall in this signal region. The number of background and signal events for a 15 tonne-year exposure for each of the benchmark scenarios and for the two xenon isotopes are shown in table 2. The upper and lower tables are for the XenonA200 and XenonB1000 benchmarks scenarios respectively. As explained in section 3, the dominant background is always from the 2νββ-decay of 136 Xe. The number of dark matter events quoted are for a dark matter particle with mass m DM = 1000 GeV and a cross-section σ 0 n = 5 × 10 −40 cm 2 . The various signal regions in the S1-log 10 (S2 b /S1) plane were taken as follows: • XenonA200, 129 Xe: 'S1 only' refers to the range 125 PE ≤ S1 ≤ 275 PE. 'S1 & log 10 (S2 b /S1)' has the additional constraint 1.15 ≤ log 10 (S2 b /S1) ≤ 1.36. Figure 11. A comparison of the discovery reach calculations using the cut-and-count approach discussed in this appendix (solid green and orange line) and the approach taken in the main body of the paper that takes into account the position of each event in the S1-log 10 (S2 b /S1) plane. The two calculations give comparable results, with the cut-and-count approach giving a limit on σ 0 n that is higher since it uses less of the information measured by our mock experiments.
These signal regions were chosen 'by eye' with reference to figure 7 so it may be possible that this calculation could be improved by choosing more optimal signal regions. As this calculation mainly serves as a cross-check of the calculation in the main body of the paper, we have not investigated this further.
In this calculation, we again employ the profile likelihood ratio described in section 4. However in this instance, the extended likelihood is We make a further simplifying assumption and consider only one normalisation of the background, rather than having an individual normalisation for each background component. In eq. (C.1), L(A BG ) is a Normal distribution that accounts for the uncertainty of the background, which we assume has an error of 4% (obtained by combining the uncertainties from the individual components). The discovery limits that follows from this calculation are shown by the solid green and orange lines in figure 11. Also show by the red and blue dashed lines are the discovery limits from the main body of the paper. As we would naively expect, the cut-and-count discovery limit on σ 0 n is higher than the result in the main body of the paper. This is because the cut-and-count calculation uses less of the information measured by our mock experiments.

JCAP05(2016)033
However, the difference is relatively small (about a factor of two) giving us confidence in the more involved discovery limit calculation presented in the main body of the paper.

D Comparing neutron-only and proton-only spin-dependent constraints
In the main body of the paper we only compared the discovery reach against the exclusion limits from two-phase xenon experiments, which are sensitive to the spin-dependent interaction with the neutron. In this appendix we discuss the limits from experiments that are sensitive to the spin-dependent interaction with the proton. A comparison of the neutron-only and proton-only spin-dependent limits are necessarily model dependent. The cross-sections from the axial-vector interaction are related by The factors g q are the coupling constants of the axial-vector mediating particle with the quarks and come from the dark matter model. When g q is the same for all quarks, σ 0 p = σ 0 n so the limits can be trivially compared. Arguably a better motivated choice is g u = −g d = −g s . This is the case with the Z-boson in the Standard Model and is therefore applicable to well-motivated dark matter candidates such as the neutralino. In this scenario, we find that σ 0 n = 0.75 × σ 0 p . The strongest direct detection constraints on the proton-only spin dependent crosssection are from PICO2L [83] and PICO-60 [84] experiments. These are plotted as the dotdashed green and purple lines in figure 12. In this figure, we have assumed that σ 0 n = 0.75×σ 0 p and rescaled the published PICO limits appropriately. By comparing with the LUX limit (grey long-dashed line) we observe that the LUX limit is always below the PICO limits. This justifies our choice of only showing the LUX limits in the main body of the paper since this experiment is the most constraining direct detection experiment for this type of interaction. There are also constraints on the proton-only spin-dependent cross-section from the IceCube search for neutrinos originating from the Sun [85]. These limits are valid when the dark matter capture and annihilation rates are in equilibrium (which must be checked on a model-by-model basis). The limits also depend on the dark matter annihilation final-state. The strongest limits are for annihilation to W + W − (brown-dashed line in figure 12), weaker limits are obtained for annihilation tobb (orange-dashed line in figure 12), while no limits exist for annihilation to first and second generation quarks and leptons. As with the PICO limits, we have assumed that σ 0 n = 0.75 × σ 0 p and rescaled the published limits appropriately. In this instance, we see that the IceCube W + W − limit is comparable to the discovery reach from a 15 tonne-year two-phase xenon experiment for a dark matter mass below ∼ 1000 GeV. However, given the model-dependent assumptions that enter the IceCube constraints, it is possible to consider scenarios in which the IceCube limit does not apply. In comparison, the LUX limit is model independent.  Figure 12. A comparison of the proton-only spin-dependent limits from PICO-2L, PICO-60 and Ice-Cube with the exclusions limit from LUX and the discovery reach for inelastic scattering off the 129 Xe and 131 Xe isotopes. We have assumed that σ 0 n = 0.75 × σ 0 p , the relationship that holds, for instance, with the Z-boson in the Standard Model, and rescaled the published PICO and IceCube limits appropriately. The IceCube limits are comparable to the discovery reach over part of the mass range but are rather model-dependent (see text for details) and it is straightforward to envisage scenarios in which these limits do not apply. The LUX limit is model-independent and is always below the PICO limits, justifying our choice to show only the LUX limits in the main body of the paper.

E The discovery reach with a halo model from EAGLE
We assumed the Standard Halo Model (SHM) in all of the calculations in the main body of the paper. This is the simplest and canonical halo model used by the direct detection community. In this appendix we investigate how our results change when another halo model is used. It is well known that the choice of halo can result in large deviations in the scattering rate for scattering processes that involve an inelastic transition. This is because these processes typically probe dark matter particles that are towards the tail of the velocity distribution [86,87]. The alternative halo that we consider in this appendix is from the EAGLE simulation [88,89], which is a state-of-the-art simulation of galaxies that contains the effects of dark matter and baryons. The simulation contains a vast number of halos so various criteria must be used to select those that are Milky Way-like. This procedure was implemented in [90,91]. The quantity which enters the event rate is the velocity integral Here we use the velocity integral that deviates most from the SHM while passing all of the selection criteria from the EAGLE HR simulation run in [91]. This is shown as the blue solid line in figure 13, where it is compared with the halo integral from the Standard Halo Model, the dashed black line. The main difference is that the EAGLE halo contains many more particles at higher values of v min . In figure 14 we show the discovery reach when using the EAGLE halo model and compared with the discovery reach from the SHM. As we would naively expect based on [86,87], we find that the largest difference between the halo models is at low mass and when the inelastic splitting is large. This is simply because the EAGLE halo has more particles with  Figure 14. A comparison of the discovery reach assuming the Standard Halo Model (dashed red and blue lines) and the EAGLE halo model (solid green and orange lines) for inelastically scattering off the 129 Xe and 131 Xe isotopes. The left and right panels show the results for the XenonA200 and XenonB1000 benchmark scenarios assuming a 15 tonne-year exposure. Also shown in the LUX exclusion limit (black dashed line) assuming the EAGLE halo model and the XENON1T projected limit assuming the SHM. Although the discovery reach assuming the SHM and EAGLE halo models show differences at smaller values of the dark matter mass, the conclusions drawn from this figure are essentially the same as in figure 8, showing that our overall results are not particularly sensitive to the choice of halo model. a higher speed, so an experiment observes more events when the dark matter mass is below approximately 100 GeV. We have also recalculated the LUX exclusion limit from their 2013 search [45] for dark matter that elastically scatters with the xenon isotopes with the EAGLE halo (we were not able to recalculate the XENON1T projected limit since we do not know all of the assumptions entering its calculation). The calculation of this limit is described in [92]. Overall, we see that our conclusions remain essentially unchanged; the discovery reach still lies between the LUX limit and XENON1T projected limit. The small difference is the that the SHM results in a slightly more conservative discovery reach, with the EAGLE discovery reach extending to slightly smaller values of the dark matter mass.