Search for inelastic scattering of WIMP dark matter in XENON1T

We report the results of a search for the inelastic scattering of weakly interacting massive particles (WIMPs) in the XENON1T dark matter experiment. Scattering off $^{129}$Xe is the most sensitive probe of inelastic WIMP interactions, with a signature of a 39.6 keV de-excitation photon detected simultaneously with the nuclear recoil. Using an exposure of 0.89 tonne-years, we find no evidence of inelastic WIMP scattering with a significance of more than 2$\sigma$. A profile-likelihood ratio analysis is used to set upper limits on the cross-section of WIMP-nucleus interactions. We exclude new parameter space for WIMPs heavier than 100 GeV/c${}^2$, with the strongest upper limit of $3.3 \times 10^{-39}$ cm${}^2$ for 130 GeV/c${}^2$ WIMPs at 90\% confidence level.

of inelastic WIMP scattering with a significance of more than 2σ.A profile-likelihood ratio analysis is used to set upper limits on the cross-section of WIMP-nucleus interactions.We exclude new parameter space for WIMPs heavier than 100 GeV/c 2 , with the strongest upper limit of 3.3 × 10 −39 cm 2 for 130 GeV/c 2 WIMPs at 90% confidence level.

I. INTRODUCTION
A wealth of astrophysical and cosmological evidence points towards the existence of dark matter [1,2].Of the many postulated candidates for dark matter, the weakly interacting massive particle (WIMP) is particularly wellmotivated, and would be expected to have directly detectable interactions with baryonic matter [3,4].A wide variety of experiments have searched for such an interaction, but a convincing signal is yet to be observed [5][6][7].
In this work, we present a search for the inelastic scattering of WIMPs off nuclei in the XENON1T experiment [8].
The main focus of direct detection searches is usually the elastic scattering of WIMP dark matter off target nuclei, aiming to detect the O(10 keV) recoiling nucleus [9][10][11][12].This work concentrates instead on inelastic scattering, where the target nucleus is left in an excited state after the interaction.The subsequent de-excitation of the nucleus produces a characteristic γ-ray which is detected as an additional energy deposit in the detector, alongside the kinetic energy transferred to the nucleus.The inelastic scattering of dark matter discussed here is not to be confused with so-called inelastic dark matter models, in which the dark matter particle itself can be excited during interactions [13].
We concentrate on scattering off a particular isotope: 129 Xe.It has the lowest-lying excited state of any xenon isotope as well as a relatively high natural abundance of 26%, meaning that it dominates the expected rate of inelastic WIMP scattering in xenon-based detectors.The first excited 3/2 + state lies 39.6 keV above the 1/2 + ground state and has a half-life of 0.97 ns [14].The structure functions for inelastic scattering off 129 Xe were calculated in [15], and the corresponding recoil spectra are used for this work.An additional channel exists in 131 Xe, with a slightly lower isotopic abundance of 21%.However, because the first excited 1/2 + state lies 80.2 keV above the 3/2 + ground state (half-life of 0.48 ns) [16], this channel is suppressed and not considered in this anal-ysis [17].
Direct detection experiments such as XENON1T tend to focus on searches for elastic WIMP scattering, where they have the highest sensitivity to most interaction models.This is due to their ability to discriminate nuclear recoil (NR) WIMP signals from electronic recoils (ERs), which constitute the majority of the backgrounds.As ERs constitute the majority of events in the energy range relevant for this search, the γ-ray component of the signal for inelastic scattering worsens this background discrimination.That said, inelastic scattering is not only a complementary search strategy.An observation of WIMP inelastic scattering would also constrain the properties of the WIMP beyond what is possible with elastic scattering only.In addition, the inelastic signal would only be detectable for spin-dependent interactions, whereas an elastic scattering signal could be detected for both spin-dependent and spin-independent interactions [17].Finally, there are some interaction models for which the inelastic channel is more sensitive [18].

II. THE XENON1T DETECTOR
Located in the Gran Sasso National Laboratory (LNGS) in central Italy, the XENON1T detector was a cylindrical dual-phase (liquid/gas) xenon time projection chamber (TPC) [8], operated between 2016 and 2018.The active target was 96 cm in diameter and 97 cm in length, enclosing 2 t of ultra-pure xenon.This was viewed from above and below by two arrays of 127 and 121 Hamamatsu R11410-21 photomultipler tubes (PMTs), respectively [19,20].An additional 1.2 t of xenon lay between the TPC and the inner cryostat vessel, providing a layer of passive shielding.Three electrodes, the cathode located near the bottom of the TPC, the gate just below the liquid surface, and the anode just above it, established the electric fields necessary for the detector operation.
The TPC itself was housed in a vacuum-insulated cryostat, which in turn was housed in a large tank ∼ 10 m in diameter and height.This tank was filled with deionised water and instrumented with 86 Hamamatsu R5912 PMTs as an active water-Cherenkov muon veto [21].A service building adjacent to the water tank holds supporting infrastructure elements such as the cryogenic system, xenon purification, distillation [22], and storage, data acquisition [23], and slow control.The materials forming the cryostats and all detector components were selected after a rigorous material screening campaign to ensure high radiopurity [24,25].
When a particle interacts in the liquid target, scintillation light is produced and xenon atoms are ionized.The scintillation light is detected by the PMTs and forms the S1 signal.The drift field between the cathode and gate drifts the liberated electrons up towards the liquid surface, where the much stronger extraction field between the gate and anode extracts the electrons into the gas and causes energetic collisions with the xenon atoms which results in further scintillation.This second, amplified signal constitutes the ionisation signal or S2, which is proportional to the number of extracted electrons.The full reconstruction of the interaction vertex is achieved through the combination of the time between S1 and S2 signals (providing the z coordinate) and the illumination pattern on the top PMT array (providing the (x, y) coordinates).The energy deposited in the interaction is reconstructed as a linear combination of the scintillation and ionisation signals [26].

III. DATA ANALYSIS
An inelastic scattering event's signature comes from two distinct energy depositions: the NR itself and a subsequent ER from the 39.6 keV γ-ray produced by nuclear de-excitation.We assume that the two energy depositions are effectively simultaneous, since the 0.97 ns half-life of the excited state is much shorter than the 10 ns time resolution of the data acquisition.The mean photon-absorption length at this energy in liquid xenon is 150 µm [27].This is much longer than the electron-ion recombination length scale of 4.6 µm [28].We therefore assume that the liquid xenon response to each deposition is unaffected by the other and that they can therefore be treated as independent.We also note that a 150 µm difference in the depth of the two interactions corresponds to a difference in arrival time for the two S2s on the order of 0.1 µs.As this is substantially smaller than the width of a typical S2 signal (O(1 µs)), we expect inelastic scattering events' S2s to appear like those from a single energy deposition.

A. Event selection
Data selection follows [9] and [29] very closely, with extended details and discussion in [30].Data were selected from the period between 2 February 2017 and 8 February 2018, when the detector conditions were stable and the drift field was constant.Because calibration signals from 83m Kr are very close to the signal region of interest in this search, data taken shortly after calibration periods using this isotope were removed, resulting in a total of 249.07 days of livetime.
From this initial data selection, additional quality requirements are imposed.To limit possible bias, all events near the signal region were blinded for this work prior to the start of this analysis.The same fiducial volume was used as in [9], containing 1316 kg of xenon.Inside this volume, events are selected that are consistent with a well-reconstructed single-site interaction and do not occur shortly after a high-energy event, where the delayed extraction of electrons [31] produces additional S2 signals that can adversely affect the event reconstruction.Events are required to have an energy deposition between 15 and 71 keV, which forms the region of interest for this analysis, based on the expected signal and backgrounds (see sections III B and III C).The combined efficiency of the selection requirements to inelastic WIMP signals is estimated at (93.6 ± 0.7)%, determined using either simulated events or calibration data.
Correction factors are applied to the measured event quantities to account for systematic variations throughout the detector volume [30].The most significant are corrections for the spatially-varying light collection efficiencies for the ionisation and scintillation signals, and the attenuation of the ionisation signal from charge loss due to electron attachment onto electronegative impurities in the liquid xenon.An additional correction accounts for inhomogeneity in the electric fields.The corrections transform the S1 signal into cS1, and the S2 into cS2.As the S2 signal is generated only a few cm below the top PMT array, the response of this set of photosensors is significantly non-uniform.Thus, only the S2 light seen by the bottom PMT array (cS2 b ) is used for energy reconstruction.The analysis is performed with the corrected quantities cS1 and cS2 b .

B. Expected signal
The response of XENON1T to an inelasticallyscattering WIMP is modelled using a combined datadriven and Monte-Carlo-driven procedure.The two energy depositions -the NR from the collision and ER from the subsequent gamma ray -are considered independently.The complete response of XENON1T to an inelastic scatter is the sum of these two models.
The response to the NR is modelled analogously to previous elastic scattering analyses of XENON1T data.The mass-dependent WIMP spectra of [15] are used together with the XENON1T response model described in [32], based on a fit to calibration events induced by neutron sources.Since the inelastic scatter is detected as a single event, we do not apply detector efficiencies to the separate signal components, but to the final combined signal model.
The ER signal is assumed to be well described by a two-dimensional (2D) Gaussian function in the plane of cS1 and cS2 b [26].The total (ER plus NR) response model for inelastic neutron scattering is fit to calibration data from a D-D neutron generator [33], using a binned likelihood maximisation procedure, to extract the five parameters needed to describe this Gaussian.The NR response is fixed for the fit, and determined using the method described above for WIMPs.A background of pure ER events is present in this calibration data, predominantly due to decays of 214 Pb.We therefore fit an  additional component, whose shape is assumed to be the same as the ER band in 220 Rn calibration data (see section III C).The best fit model for inelastic neutron scattering is shown in figure 1 along with the neutron generator calibration data used for the fit.The goodness-of-fit p-value is 0.22, computed from the likelihood compared to its distribution for toy data.
The ER model produced with these best-fit parameters is then combined with the NR model for a particular WIMP mass.After considering the selection efficiencies, this gives the total expected response to such a WIMP inelastically scattering in XENON1T.

C. Backgrounds
Four primary backgrounds contribute to the signal region of this analysis.These are decays of 214 Pb (a daughter of 222 Rn), residual contamination of the calibration isotope 83m Kr, and peaks from 124 Xe and 125 I [29].Additional contributions from elastic scattering of solar neutrinos off electrons, decays of 85 Kr, 2ν2β decays of 136 Xe, and Compton scatters of gammas emitted by detector components are suppressed by at least one order of magnitude and not considered.
The backgrounds presented here are modelled in the parameter space described by the ratio cS2 b /cS1 and energy.Although ER and NR events have different energy scales, we use the ER-equivalent energy.The 1σ and 2σ contours of each model are shown along with the expected signal for a 100 GeV/c 2 WIMP in figure 2.

214 Pb
The isotope 222 Rn continuously emanates from all detector materials due to trace contamination of 238 U [25].Itself a noble element, radon is unaffected by purification techniques which target electronegative impurities such as oxygen and water.Its half-life of 3.8 days is sufficiently long for it to mix uniformly throughout the active volume [34], where it decays first to 218 Po (halflife 3.1 min) then 214 Pb (half-life 27 min).Of the decays of 214 Pb, 11% are directly to the ground state of 214 Bi and have no accompanying gamma emission [35].These "naked" beta decays can therefore deposit any energy up to the Q-value of 1.02 MeV, and some fraction of these fall within the energy region of interest in this analysis.
The spectral shape of this background is modelled using Monte Carlo simulations.The cS2 b /cS1 distribution is based on 220 Rn calibration data with decays of 212 Pb.

83m Kr
A low rate of 83m Kr -on average 10 −4 times the rate during calibrations -persisted in the detector even considerably after the injection of this isotope for calibration purposes, possibly due to trace contamination from its parent isotope 83 Rb in the xenon purification loop.Events where the two constituent decays are unresolvable result in a single signal at 41 keV, which lies in the energy region of interest.We use 83m Kr calibration data, after applying the analysis selection criteria, to model the distribution of this background.The ER background stemming from 124 Xe and 125 I is modelled in two ways.Each of these two isotopes provides two peaks in the region of interest, depending on which electron shells are involved in the decay. 124Xe, decaying via two-neutrino double electron capture (2νECEC), produces peaks at 64.3 keV (36.7 keV) with branching ratios ∼ 76% (∼ 23%) via capture from the KK (KL) electron shells (here denoted 2νECEC-KK and 2νECEC-KL, respectively) [36]. 125I is produced by neutron capture on 124 Xe primarily during neutron calibrations, and decays via electron capture (EC).This produces peaks at 67.3 keV (40.4 keV) with branching ratios 80.1% (15.6%) via capture from the K (L) shell (here denoted EC-K and EC-L, respectively) [37].
The models for these backgrounds are created in a different manner than for the previous two, since no directly comparable calibration data is available.The 2νECEC-KK and EC-K events lie relatively far from the signal region, so a model consisting of one-dimensional Gaussian peaks in energy is sufficient.
For the 2νECEC-KL and EC-L decays, a onedimensional model would have a significant adverse impact on the sensitivity.In the following we describe the production of the two-dimensional model for 2νECEC-KL; the model for EC-L is produced analogously.The two-dimensional model is formed by scaling the light and charge yields of the 83m Kr calibration peak.The contributions from X-rays and Auger electrons to the total energy deposition of 2νECEC-KL decays are determined using RELAX [38].The NEST package [39] is then used to simulate the light and charge yields of liquid xenon to pure betas and pure gammas at the total energy deposition.These are linearly combined based on their fractional contributions to estimate the total light and charge yield for 2νECEC-KL.The predicted light and charge yields for 83m Kr decays are directly obtained from NEST.To estimate the cS1 (cS2 b ) distribution for 2νECEC-KL, the distribution observed in 83m Kr calibration data is scaled by the ratio between the 2νECEC-KL and 83m Kr light (charge) yields.

D. Statistical interpretation
For each WIMP mass considered between 20 GeV/c 2 and 10 TeV/c 2 , a binned profiled likelihood is used to constrain the cross-section of inelastic WIMP interactions.The binning structure used to evaluate the likelihood, shown in figure 2, is chosen in order to optimise sensitivity while preserving asymptotic properties of the likelihood [40].This explains the larger bins used in the tail regions of the background, at high and low values of cS2 b /cS1.Specifically, we ensure that a minimum of five background events are expected in every bin.At energies above 55 keV, where 2νECEC-KK and EC-K backgrounds are present, bins are only a function of energy, since these peaks are not modelled in two dimensions.Events between 15 and 29 keV are taken together to constrain the 222 Rn rate.This is taken into account using as a single ancillary likelihood term that is combined with the binned likelihood.
The data-derived backgrounds from 222 Rn and 83m Kr have an uncertainty on the fraction of events expected in each bin due to the finite statistics of the calibration data used to create the model.These statistical uncertainties can be treated with a simultaneous fit to the calibration and the science data.We follow the method described in [41] which makes it possible to incorporate this into the fitting procedure in a computationally efficient way.
During the operation of XENON1T, a very high number of events from the decay of 83m Kr were recorded during the regular calibrations.Therefore, in every bin, the statistical uncertainty arising from 83m Kr is never more than 0.4% of that from the 222 Rn model and can thus be neglected.The effect of statistical fluctuations in the Monte-Carlo data used to produce the 222 Rn model is also ignored -at any given energy the statistics is between 700 and 800 times greater than the 220 Rn calibration data.This means that the simultaneous fit is only to one calibration source ( 220 Rn) in addition to the science data.In this case the method of [41] becomes analytic, and no additional degrees of freedom must be introduced into the minimisation routine compared to fitting the science data alone.
Uncertainties on the mean energy and resolution of the 2νECEC-KK and EC-K peaks are incorporated as independent nuisance parameters in the fit, each with a Gaussian constraint term.The uncertainty on the charge yield for the 2νECEC-KL and EC-L peaks is also included with a Gaussian constraint.To constrain these parameters, the charge yield in the ER band at the energy of the 83m Kr peak is compared to the yield predicted by scaling the 83m Kr peak itself using the same procedure as for producing the models, described above.The difference between predicted and measured charge yields in this case is taken as an uncertainty on the predicted yields for the 2νECEC-KL and EC-L peaks.Because the physical processes behind each peak are very similar, we assume that the uncertainties are correlated and a single nuisance parameter is used to vary the charge and light yields for both peaks simultaneously.
The effect on the signal model of eighteen parameters is considered.These are the five parameters needed to describe the 2D Gaussian part of the model and the thirteen parameters from [32] which are relevant to the NR model.The importance of each of these is assessed by determining how much the sensitivity changes when that parameter is varied by ±1σ with the remaining parameters fixed to their nominal values.The mean cS2 b of the ER part of the signal is found to be the most important, especially at low WIMP masses (3.4% effect at 20 GeV/c 2 ), and this parameter is included in the likelihood fit with a Gaussian constraint.The others are combined into a single uncertainty on the signal rate, The inset shows the same data rebinned according to the ratio between the expected number of signal events from a 50 GeV/c 2 WIMP and the expected number of background events, from the background-only best fit.calculated as the sum in quadrature of their effect on the sensitivity.This is also combined with the uncertainty on the efficiency of the selection criteria to form a single nuisance parameter which we term the "effective efficiency".A summary of the constrained parameters is shown in Table I.

IV. RESULTS AND DISCUSSION
The distribution of events observed in the 0.89 tonneyear exposure after unblinding is shown in figure 3. A total of 7392 events are observed.The number of events expected from each source of background, for the best fit models, and the best fit number of signal events, is detailed in table II.Based on the profiled likelihood analysis, no significant evidence is found for spin-dependent inelastic WIMP-nucleon scattering in the search data.A pre-unblinding decision was made to report only an upper limit if the discovery significance p-value was less than 0.003 (3σ).The highest significance is for 50 GeV/c 2 WIMPs, where the p-value for the background-only hypothesis is 0.023, from the log likelihood ratio.A 90% confidence level upper limit is therefore placed on the cross section.Under the background-only scenario, we obtain a goodness-of-fit p-value of 0.055, by comparing the χ 2 goodness of fit to its distribution for simulated data.
The result is shown in figure 4, along with the expected 1σ and 2σ range of upper limits.We set the strongest upper limit for WIMP masses greater than 100 GeV/c 2 , reaching 3.3×10 −39 cm 2 at a WIMP mass of 130 GeV/c 2 .The sensitivity to lower-mass WIMPs (up to about 30 GeV/c 2 ) is severely affected due to the similarity between their signal shape and the 83m Kr background.For the lowest masses considered, this effect weakens the sensitivity by approximately a factor of three; above 130 GeV/c 2 the effect is approximately 10%.Our limit is weaker than the upper 1σ quantile of the expected range of limits for all WIMP masses.This is due to a slight over-fluctuation in the bins with the most sensitivity to WIMPS, as seen in the inset of figure 3. We also com-  The expected range at 1σ and 2σ is shown in green and yellow, respectively.Our result is compared to upper limits obtained by the XMASS collaboration [42] (blue dot-dashed line) and XENON100 [43] (orange dashed line).
pare the limit to that reported by the XMASS collaboration in [42], obtained from a 0.72 tonne-year exposure of XMASS-I.
These results complement previously reported searches for elastic WIMP scattering [9,44].The sensitivity is lower for most WIMP-scattering scenarios.However, a positive detection of inelastic scattering would place stronger constraints on the properties of the dark matter than a detection of elastic scattering alone.
Next-generation experiments including XENONnT, which is currently being commissioned, are expected to have an ER background rate around five times lower [45].In addition, a much larger sensitive region will mean that larger exposures can quickly be obtained.This will allow even smaller cross-sections of inelastic scattering to be probed.

FIG. 1 .
FIG. 1. Model for inelastic scattering of neutrons (pink 1σ, 2σ and 3σ contours) compared to calibration data with the neutron generator deployed.The density of calibration data is represented by the shading inside the 2σ contour, while outside individual events in the calibration data are shown as green dots.The white unshaded region indicates the region in which the fit is performed.Energy contours are shown as grey lines.

FIG. 2 .
FIG. 2. Region interest showing the four nominal background models: 214 Pb in grey, 83m Kr in green, 124 Xe in pink and 125 I in blue; and the expected signal for a 100 GeV/c 2 WIMP in red.In each case the darker and lighter shading shows the 1σ and 2σ region, respectively.The cS2 b /cS1 distribution for 214 Pb is used to visualise the 2νECEC-KK and EC-K peaks, which are only modelled in energy.The dotted lines indicate the boundaries of the bins used to perform the likelihood fit.

FIG. 3 .
FIG.3.Observed events in the 0.89 t exposure, shown as error bars.Between 29 and 55 keV, the six histograms correspond to the six cS2 b /cS1 bins as shown in figure2.Above 55 keV only a single cS2 b /cS1 bin is shown, as used for the analysis.The single bin from 15 to 29 keV is not shown, since it is instead used to constrain the background rate in the rest of the analysis region.The shaded areas show the estimated background in each bin from 214 Pb (grey), 83m Kr (green), 124 Xe (pink) and 125 I (blue), based on a background-only best fit to the science data.The solid (dotted) red line shows the signal expected from a 50 GeV/c 2 (100 GeV/c 2 ) WIMP scattering inelastically, with a cross-section of 10 −37 cm 2 , on top of the background estimation.The inset shows the same data rebinned according to the ratio between the expected number of signal events from a 50 GeV/c 2 WIMP and the expected number of background events, from the background-only best fit.

FIG. 4 .
FIG.4.Upper limit, at 90% CL, on the cross-section for inelastic scattering of WIMPs off xenon nuclei.The expected range at 1σ and 2σ is shown in green and yellow, respectively.Our result is compared to upper limits obtained by the XMASS collaboration[42] (blue dot-dashed line) and XENON100[43] (orange dashed line).

TABLE I .
Overview of parameters with a Gaussian constraint in the likelihood.The effective efficiency has a massdependent uncertainty, varying between 0.936 ± 0.033 for 30 GeV/c 2 WIMPs and 0.936 ± 0.014 for 10 TeV/c 2 WIMPs.

TABLE II .
Best fit expectation values for the number of events from 100 GeV/c 2 WIMPs and each background, under the likelihood's best fit scenario."All bins" refers to the entire region of interest (15-71 keV, 0 ≤ cS2 b /cS1 ≤ 50), and "signal region" refers to the bins containing 95% of the signal model.