XENON1T Dark Matter Data Analysis: Signal & Background Models, and Statistical Inference

TL;DR

This paper details the XENON1T experiment's detector response, background and signal models, and statistical methods, leading to the most stringent limits on WIMP-nucleon interactions for certain dark matter masses.

Contribution

It provides a comprehensive description of the detector response, background, and signal models, along with statistical inference procedures, advancing dark matter search techniques with XENON1T data.

Findings

Achieved the lowest background rate in liquid xenon detectors

Set the best limit on WIMP-nucleon cross-section for masses above 6 GeV/c^2

Developed detailed models and inference methods for dark matter detection

Abstract

The XENON1T experiment searches for dark matter particles through their scattering off xenon atoms in a 2 tonne liquid xenon target. The detector is a dual-phase time projection chamber, which measures simultaneously the scintillation and ionization signals produced by interactions in target volume, to reconstruct energy and position, as well as the type of the interaction. The background rate in the central volume of XENON1T detector is the lowest achieved so far with a liquid xenon-based direct detection experiment. In this work we describe the response model of the detector, the background and signal models, and the statistical inference procedures used in the dark matter searches with a 1 tonne$\times$year exposure of XENON1T data, that leaded to the best limit to date on WIMP-nucleon spin-independent elastic scatter cross-section for WIMP masses above 6 GeV/c$^2$.

Figures

Tables

Table 1. Table 1: Goodness of fit p-values of data-model matching on cS2 b spectra in different cS1 ranges for 220 Rn, 241 AmBe and D-D neutron generator calibration data in SR1.

data	cS1 range (PE)
	0-10	10-20	20-30	30-40	40-50	50-70
220Rn overall	0.18	0.44	0.51	0.96	0.29	0.49
220Rn reference	0.17	0.28	0.30	0.85	0.21	0.48
241AmBe	0.12	0.50	0.88	0.89	0.62	0.70
D-D	0.10	0.50	0.89	0.73	0.11	0.27

Table 2. Table 2: Main nuisance parameters in the detector response model of XENON1T. Lower and upper errors of the posterior are calculated as the difference between the median and the 15.4% and 84.6% percentiles of the posterior.

Par.	Prior		Posterior		Reference and note
	SR0	SR1	SR0	SR1
$W$	13.7$\pm$0.2		13.8$\pm$0.2		In unit of eV; Eq. 1
$⟨N_{\mathrm{ex}}/N_i⟩$	0.06 - 0.20		${0.15}_{-0.06}^{+0.04}$		Eq. 2
$g_1$	0.142$\pm$0.002		0.142$\pm$0.005		In unit of PE/photon; Eq. 10
$g_2$	11.4$\pm$0.2		11.4$\pm$0.2		In unit of PE/e-; Eq. 10, for cS2b
$p_{\mathrm{dpe}}$	0.18 - 0.24		${0.219}_{-0.023}^{+0.015}$		Eq. 10, Faham et al. (2015); Paredes et al. (2018)
1 - ${\tau }_e$/${\tau }_e^{\prime }$	0$\pm$0.04	0$\pm$0.02	0.01$\pm$0.03	0.01$\pm$0.01	Eq. 12 and 16
${\gamma }_{\mathrm{er}}$	free		0.124$\pm$0.003		Eq. 8
${\omega }_{\mathrm{er}}$	free		31$\pm$4		In unit of keV; Eq. 8
${\delta }_{\mathrm{er}}$	free		0.24$\pm$0.06		Eq. 8
$q_0$	free		${1.13}_{-0.32}^{+0.24}$		In unit of keV; Eq. 7
$q_1$	free		${0.47}_{-0.15}^{+0.18}$		In unit of keV; Eq. 7
$q_2$	free	free	0.041$\pm$0.006	0.034$\pm$0.003	Eq. 6
$q_3$	free		1.7${}_{-1.1}{}^{}{}_{}{}^{+1.3}$		In unit of keV; Eq. 6
$\alpha$	1.240$\pm$0.079		1.280$\pm$0.063		NEST parameters Lenardo et al. (2015)
$\zeta$	0.047$\pm$0.009		0.045${}_{-0.008}{}^{}{}_{}{}^{+0.009}$		NEST parameters Lenardo et al. (2015)
$\beta$	239$\pm$28		273${}_{-20}{}^{}{}_{}{}^{+24}$		NEST parameters Lenardo et al. (2015)
$\gamma$	0.0139$\pm$0.007		0.0141$\pm$0.006		NEST parameters Lenardo et al. (2015)
$\delta$	0.062$\pm$0.006		0.061$\pm$0.006		NEST parameters Lenardo et al. (2015)
$\kappa$	0.139$\pm$0.003		0.138$\pm$0.003		NEST parameters Lenardo et al. (2015)
$\eta$	3.3$\pm$0.7		3.3${}_{-0.6}{}^{}{}_{}{}^{+0.7}$		NEST parameters Lenardo et al. (2015)
$\lambda$	1.14$\pm$0.45		1.15${}_{-0.27}{}^{}{}_{}{}^{+0.35}$		NEST parameters Lenardo et al. (2015)
${ϵ}_{\mathrm{ext}}$	96%	96%	-	-	Fixed; Eq. 12
$\Delta G/G$	0.24	0.25	-	-	Fixed; Eq. 13
${\nu }_d$	0.144	0.134	-	-	In unit of cm/$\mu$s; fixed; Eq. 12
$F$	120	81	-	-	In unit of V/cm;

Table 3. Table 3: Table of parameters of the full XENON1T likelihood, best-fit values, and 1 ​ σ 1 𝜎 1\,\sigma confidence interval for a 200 GeV / c 2 200~{}\mathrm{GeV}/\mathrm{c}{}^{2} WIMP. In addition, the data-set that constrains the parameter is noted. Parameters are arranged according to whether they affect background rates, the shape of the background distributions, or the signal distribution. For the shape parameters, the models affected by the parameter are also noted. Due to the mismodelling parameter, the signal distribution, which changes with the WIMP mass, may also affect the ER model best-fit. Expectation values are written as events in each dataset; SR0 (SR1) expectation values refer to a 0.12 ​ t × y 0.12 t y 0.12~{}\mathrm{t\times y} ( 0.88 ​ t × y 0.88 t y 0.88~{}\mathrm{t\times y} ) exposure. The radiogenic neutron event rate, shared between the science runs, refers to the 1.0 ​ t × y 1.0 t y 1.0~{}\mathrm{t\times y} exposure.

Rate Parameter		Constraint	Expectation Value
SR0 AC events		Ancillary measurement	${0.051}_{-0}^{+0.035}$
SR1 AC events		Ancillary measurement	${0.42}_{-0}^{+0.29}$
SR0 220Rn  AC events		Ancillary measurement	${0.73}_{-0}^{+0.51}$
SR1 220Rn  AC events		Ancillary measurement	${2.8}_{-0}^{+1.9}$
SR0 CE$\nu$NS events		Ancillary measurement	${0.0040}_{-0.0014}^{+0.0014}$
SR1 CE$\nu$NS events		Ancillary measurement	${0.050}_{-0.017}^{+0.017}$
SR0 science data ER		SR0  data-set	${73.4}_{-8.4}^{+8.8}$
SR1 science data ER		SR1  data-set	${554}_{-24}^{+24}$
SR0 calibration data ER		SR0  220Rn  data-set	${689}_{-26}^{+27}$
SR1 calibration data ER		SR1  220Rn  data-set	${5264}_{-72}^{+73}$
Radiogenic events in SR0+SR1		Ancillary measurement	${1.44}_{-0.66}^{+0.65}$
SR0  surface events		SR0 science data	${12.5}_{-3.2}^{+3.9}$
SR1  surface events		SR1 science data	${93.7}_{-9.5}^{+10}$
Shape Parameter	Model Affected	Constraint	Value
ER photon yield	All ER models	220Rn calibration	$-{0.043}_{-0.052}^{+0.053}$
SR0 ER recombination fluctuation	SR0 220Rn, science models	SR0 220Rn data-set	${0.64}_{-0.70}^{+0.58}$
SR1 ER recombination fluctuation	SR1 220Rn, science models	SR1 220Rn data-set	${0.32}_{-0.39}^{+0.41}$
SR0 safeguard as fraction of ER	SR0 220Rn, science models	SR0 220Rn data-set	$-{0.0060}_{-0.0036}^{+0.0049}$
SR1 safeguard as fraction of ER	SR1 220Rn, science models	SR1 220Rn data-set	$-{0.0049}_{-0.0045}^{+0.0014}$
SR0 surface shape parameter	SR0 surface model	SR0 science data	${1.00}_{-0.67}^{+0.0}$
SR1 surface shape parameter	SR0 surface model	SR1 science data	$-{0.39}_{-0.22}^{+0.24}$
Signal Parameter		Constraint	Value
SR0  signal efficiency		NR model uncertainty	${1.000}_{-0.036}^{+0.036}$
SR1  signal efficiency		NR model uncertainty	${1.000}_{-0.021}^{+0.021}$
WIMP cross-section [${10}^{-45}{\mathrm{cm}}^2$]			${0.0421}_{-0}^{+0.1165}$
WIMP mass [GeV/c2]	Signal, ER mismodelling term	Fixed in analysis	200

Equations

$$N_q\sim \mathrm{Binom}\left(\epsilon /W,L\right),$$

$$\large\begin{array}[]{rcl}N_{i}&\sim&\mathrm{Binom}\left(N_{q},\frac{1}{1+\langle N_{\mathrm{ex}}/N_{i}\rangle}\right),\\ N_{\mathrm{e}x}&=&N_{q}-N_{i}.\end{array}$$

$$\large\begin{array}[]{rcl}N_{e}&\sim&\mathrm{Binom}\left(N_{i},1-r\right),\\ N_{\gamma}&=&N_{i}-N_{e}+N_{\mathrm{ex}},\end{array}$$

$$r\sim \mathrm{Gauss}\left(⟨r⟩,\Delta r\right).$$

$$⟨r⟩=1-\frac{\ln (1+N_i\varsigma /4)}{N_i\varsigma /4},$$

$$\Delta r=q_2(1-e^{-\epsilon /q_3}),$$

$$⟨r{⟩}_{er}=\left(1-\frac{\text{ln}(1+N_i{\varsigma }_{\mathrm{er}}/4)}{(N_i{\varsigma }_{\mathrm{er}}/4)}\right)/\left(1+e^{-(\epsilon -q_0)/q_1}\right),$$

$${\varsigma }_{\mathrm{er}}={\gamma }_{\mathrm{er}}{e}^{-\epsilon /{\omega }_{\mathrm{er}}}{F}^{-{\delta }_{\mathrm{er}}},$$

$$\large\begin{array}[]{ccl}\langle N_{\gamma}\rangle/\varepsilon&=&\frac{1}{W}\frac{\langle r\rangle+\langle N_{\mathrm{ex}}/N_{i}\rangle}{1+\langle N_{\mathrm{ex}}/N_{i}\rangle},\\ \langle N_{e}\rangle/\varepsilon&=&\frac{1}{W}\frac{1-\langle r\rangle}{1+\langle N_{\mathrm{ex}}/N_{i}\rangle}.\end{array}$$

$$\large\begin{array}[]{ccl}g^{\prime}_{1}(x,y,z)&=&(1+p_{\mathrm{dpe}})\cdot\epsilon_{L}(x,y,z)\cdot\epsilon_{\mathrm{QE}}\cdot\epsilon_{\mathrm{CE}},\\ g^{\prime}_{2}(x,y)&=&\epsilon_{\mathrm{ext}}G(x,y),\end{array}$$

$$\large\begin{array}[]{ccl}N_{hit}&\sim&\text{Binom}\left(N_{\gamma},\epsilon_{L}(x,y,z)\cdot\epsilon_{\mathrm{QE}}\cdot\epsilon_{\mathrm{CE}}\right),\\ N_{\mathrm{pe}}-N_{\mathrm{hit}}&\sim&\text{Binom}\left(N_{\mathrm{hit}},p_{\mathrm{dpe}}\right).\end{array}$$

$$N_{\mathrm{ext}}\sim \text{Binom}\left(N_e,e^{-z/({\tau }_e\cdot {\nu }_d)}{ϵ}_{\mathrm{ext}}\right),$$

$$N_{\mathrm{prop}}\sim \text{Gauss}\left(N_{\mathrm{ext}}G,\sqrt{N_{\mathrm{ext}}}\Delta G\right),$$

$$\large\begin{array}[]{ccl}\mathrm{S1}/N_{\mathrm{pe}}-1&\sim&\text{Gauss}\left(\delta_{\mathrm{s1}},\Delta\delta_{\mathrm{s1}}\right),\\ \mathrm{S2}/N_{\mathrm{prop}}-1&\sim&\text{Gauss}\left(\delta_{\mathrm{s2}},\Delta\delta_{\mathrm{s2}}\right),\end{array}$$

$${\vec{x}}_r\sim \text{Gauss}\left(\vec{x},{\sigma }_p\right),$$

$$\large\begin{array}[]{ccl}\text{cS1}&=&\text{S1}\frac{g_{1}}{g^{\prime}_{1}(x_{r},y_{r},z_{r})},\\ \text{cS2}&=&\text{S2}\frac{g_{2}}{g^{\prime}_{2}(x_{r},y_{r})}e^{z/(\tau_{e}^{\prime}\cdot\nu_{d})},\end{array}$$

$$q(\sigma )\equiv -2\cdot \log \frac{\mathcal{L}(\sigma ,\widehat{\hat{\mathbf{\theta }}})}{\mathcal{L}(\hat{\sigma },\hat{\mathbf{\theta }})},$$

$$p_{H_0}={\int }_{q_0^{\mathrm{obs}}}^{\infty }\left[f_{H_0}(q_0)\right]\mathrm{d}{q}_0.$$

$$\text{centering}\log {\mathcal{L}}_{\mathrm{total}}(\sigma ,\mathbf{\theta })=\text{@add@centering}$$

$$+\sum _{\mathrm{SR}}\log {\mathcal{L}}_{\mathrm{SR}\mathrm{cal}}(\mathbf{\theta })$$

$$+\sum _m\log {\mathcal{L}}_m({\theta }_m),$$

$${\mathcal{L}}_{\mathrm{SR}\mathrm{sci}}(\sigma ,\mathbf{\theta })=$$

$$\times \prod _{i=1}^{N_{\mathrm{SR}}}\left[\sum _c\frac{{\mu }_c(\sigma ,\mathbf{\theta })}{{\mu }_{\mathrm{tot}}(\sigma ,\mathbf{\theta })}\cdot f_c({\vec{x}}_i|\mathbf{\theta })\right],$$

$${\mu }_{\mathrm{tot}}(\sigma ,\mathbf{\theta })\equiv$$

$${\mathcal{L}}_m({\mu }_c)=\mathrm{Gaus}({\hat{\mu }}_c|{\mu }_c,{\sigma }_c).$$

$${\mu }_{\mathrm{sig}}=\sigma \cdot ϵ\cdot {\mu }_{\mathrm{ref}}(M)/{\sigma }_{\mathrm{ref}}$$

$${\mathcal{L}}_{\mathrm{m}}(ϵ)=\mathrm{Gaus}(1|ϵ,{\sigma }_{ϵ}(M)),$$

$${\mathcal{L}}_{\mathrm{m}}({\mu }_{\mathrm{AC}})=\mathrm{Uniform}({\mu }_{\mathrm{AC}}|\; 0.6\cdot {\hat{\mu }}_{\mathrm{AC}},{\hat{\mu }}_{\mathrm{AC}}).$$

$$f_{\mathrm{ER}}^{\prime }\equiv b\cdot \left[(1-a_{\mathrm{mm}})\cdot f_{\mathrm{ER}}(\vec{x}|{\mathbf{\theta }}_{\mathrm{ER}})+a_{\mathrm{mm}}\cdot f_{\mathrm{sig}}(\vec{x})\right],$$

Full text

Dark Matter, Direct Detection, Xenon, Data Analysis

I Introduction

The existence of dark matter (DM) and its making up about 26% of the mass-energy of the Universe Aghanim et al. (2018) is indicated by a wide range of astronomical and cosmological observations. Direct detection experiments, which search for DM particles interacting with ordinary matters in a terrestrial detector target, have not yet yielded unequivocal evidence for dark matter Akerib et al. (2017); Aprile et al. (2016a); Cui et al. (2017); Abea et al. (2018); Petricca et al. (2017); Agnese et al. (2015a). The XENON1T experiment Aprile et al. (2017a), located in the INFN Laboratori Nazionali del Gran Sasso, Italy, primarily searches for Weakly Interacting Massive Particles (WIMPs), which could scatter elastically off xenon atoms. Using a 1 tonne$\times$year exposure and a nuclear recoil (NR) energy range from 4.9 keV to 40.9 keV, XENON1T has set upper limits on the cross section of spin-independent elastic scattering with a minimum of 4.1$\times$10*-47* cm2 for a 30 GeV/c2 WIMP Aprile et al. (2018a). These are the most stringent constraints set on this interaction for WIMP masses above 6 GeV/c2. The XENON1T experiment has achieved the lowest background rate among liquid xenon (LXe) detectors to date.

The dual-phase time projection chamber (TPC) used by the XENON1T detector allows for the reconstruction of the deposited energy and the three-dimensional position of interactions in the active liquid xenon target. The observable signals are the scintillation (S1) and ionization (S2) signals produced by energy depositions. The longitudinal ($z$) position is reconstructed using the time difference between the prompt S1 signal and the S2 signal, which is produced by electroluminescence in gaseous xenon after electrons, drifted upwards by an electric field, get extracted from the liquid into the gas. Both signals are observed by arrays of photomultiplier tubes (PMTs) arranged at the top and bottom of the detector. The position in the ($x$,$y$) plane is reconstructed using the S2 signal pattern in the upper PMT array. Background from radioactivity in detector materials can be rejected to a large extent by selecting a three-dimensional fiducial region within the active volume. In addition, the S2-S1 ratio can be used to discriminate between NRs from WIMPs and neutrons and electronic recoils (ERs) from $\gamma$ and $\beta$, which constitute the major backgrounds of the XENON1T experiment. More details on TPC working principles and the XENON1T TPC can be found in Aprile et al. (2017a).

The XENON1T data analysis can be divided into two parts. The first part includes event reconstruction, signal corrections, and event selection, and is reported in Aprile et al. (2019a). The second part includes the detector response model, the background and WIMP signal models, and the statistical inference, and is presented in detail in this manuscript. These models and techniques are used in the XENON1T DM searches Aprile et al. (2017b, 2018a). The detector response model, which will be presented in Section II, describes how an ER or a NR energy deposition is reconstructed in the TPC. The fit of the detector response model to calibration data provides the ER and NR background models and the signal model described in Section III, which also considers background models constructed using data-driven methods. Lastly, the statistical inference is presented in Section IV. A summary is then given in Section V.

II Detector Response Model

Understanding the conversion from the deposited energy to the observed S1 and S2 signals is critical for interpreting the results of DM searches in XENON1T. The energy region of interest is in the range of a few keV to several tens of keV in searches for elastic scatters between WIMP and xenon nuclei. The conversion of deposited energy to S1 and S2 in this region is non-linear, with fluctuations due to the scintillation and ionization processes in LXe and due to the detector reconstruction.

The model of the XENON1T detector response to ERs and NRs is based on simulations which include a comprehensive description of the signal production process and detailed characterizations of detector detection and reconstruction effects. The model is constrained by a Bayesian simultaneous fit to ER and NR calibration data, which allows to use all available information and treat correlated detector uncertainties coherently.

II.1 Basic Signal Response in Liquid Xenon

The intrinsic signal response model in LXe follows the approach used in the Noble Element Simulation Technique (NEST) model Szydagis et al. (2011); Lenardo et al. (2015). There are three forms of energy deposition in LXe: thermalization of the recoiling particle, excitation of xenon atoms, and ionization of xenon atoms. The thermalization energy loss is undetectable in the XENON1T detector. The number of detectable quanta $N_{q}$ is the sum of the number of excitons $N_{\mathrm{ex}}$ and ion-electron pairs $N_{i}$, and can be used to reconstruct the deposited energy $\varepsilon$. It follows a Binomial fluctuation due to the potential energy loss to thermalization,

[TABLE]

where $L$ is the Lindhard factor expressing the fraction of energy loss to heat and $W$ (13.7$\pm$0.2 eV from a global fit Szydagis et al. (2011)) is the average energy required to create either an exciton or ion-electron pair in LXe. Negligible energy is lost to thermalization in an ER as the mass of the recoiling electron is much smaller than the xenon nucleus. In NRs, the recoiling xenon atom transfers kinetic energy through elastic scattering off surrounding xenon atoms, resulting in a Lindhard factor Lindhard et al. (1963) of 0.1-0.2 in LXe. The field- and energy-dependence of the Lindhard factor are parametrized following the NEST model Lenardo et al. (2015).

The exciton-to-ion ratio $\langle N_{\mathrm{ex}}/N_{i}\rangle$ is related to the excitation and ionization cross sections of recoiling particles on xenon atoms. For ERs, it is assumed to be constant, and is given a uniform prior ranging from 0.06 to 0.20 Szydagis et al. (2011). For NRs, it is parametrized as a function of the deposited energy and the electric field strength $F$ in the active volume, following Lenardo et al. (2015), and is in the range of 0.7 to 1.0 for NR energies from 5 to 40 keVnr under a field of 81 kV/cm. The binomial fluctuations of $N_{ex}$ and $N_{i}$ can be written as

[TABLE]

Ionization electrons have a probability $1-r$ to escape the cloud of ion-electron pairs, where $r$ is referred to as the recombination fraction,

[TABLE]

where $N_{\gamma}$ and $N_{e}$ are the number of photons generated by de-excitation of the initial excitons and ion-electron recombination, and of the escaping electrons, respectively. Due to detector effects, such as field non-uniformity, and intrinsic fluctuations Akerib et al. (2016a), the recombination fraction $r$ fluctuates, and is modeled as a Gaussian distribution,

[TABLE]

The mean recombination fraction $\langle r\rangle$ depends on the deposited energy and the electric field, and is described by the Thomas-Imel (TI) box model Thomas and Imel (1987),

[TABLE]

where $\varsigma$ is the field-dependent TI model parameter. $\Delta r$ is the recombination fluctuation, parameterized as

[TABLE]

where $q_{2}$ and $q_{3}$ are free parameters. The parametrization is empirically chosen to take into account both the fact that a constant recombination fluctuation is observed with deposit energy larger than 2 keV and the assumption that the recombination fluctuation diminishes as deposit energy goes to zero.

For NR, the TI box model, together with the Lindhard factor, has been shown to match data well Lenardo et al. (2015). The recombination fluctuation for NRs was shown to be much smaller than the statistical fluctuations on recombination Aprile et al. (2013) induced by Eq. 2 and 3, and is set to 0 in the detector response model. For the parameterization of the NR recombination fraction $\langle r\rangle_{\mathrm{nr}}$, the TI parameter $\varsigma_{nr}$ is expressed, following Lenardo et al. (2015), as a power law function for the field dependence. For low energy ERs (roughly above 3 keV and below 10 keV), recent measurements Goetzke et al. (2017); Akerib et al. (2016a); Huang (2015); Boulton et al. (2017) indicate that the TI box model cannot fully describe the recombination process. Therefore, we use a modified TI box model for the ER recombination fraction $\langle r\rangle_{\mathrm{er}}$,

[TABLE]

[TABLE]

where the Fermi-Dirac term in Eq. 7 and the exponential term in Eq. 8 were empirically added to the TI box model to account for the deviation of measurements in the $<$3 keV and $>$10 keV energy ranges, respectively. Similarly to NR, the field dependence of the ER TI box parameter $\varsigma_{\mathrm{er}}$ follows a power law as introduced in NEST Lenardo et al. (2015). The free parameters $q_{0}$, $q_{1}$, $\gamma_{\mathrm{er}}$, $\omega_{\mathrm{er}}$, and $\delta_{\mathrm{er}}$ are obtained by matching the detector response model to XENON1T data, without any additional constraints. It is worth noting that NEST model has been updated based on a global fit using recent measurements Szydagis et al. (2018), and is compatible with this work in the energy region of interest. Fig. 1 shows the mean photon and charge yields as a function of energy for NR and ER, respectively, together with the measurements from Aprile et al. (2005, 2006, 2009, 2013); Plante et al. (2011); Sorensen et al. (2009); Manzur et al. (2010); Akerib et al. (2016b) for NR and Aprile et al. (2018b); Huang (2015); Akerib et al. (2016a); Boulton et al. (2017) for ER. The mean photon $\langle N_{\gamma}\rangle/\varepsilon$ and charge yields $\langle N_{e}\rangle/\varepsilon$ are defined as

[TABLE]

The calibration of low energy ERs in XENON1T is performed using an internal 220Rn source Aprile et al. (2017c). The energy spectrum of $\beta$-decays from 212Pb, one of the progenies of 220Rn, is similar to the dominant ER background, from $\beta$-decays of 214Pb originating from 222Rn emanation, in the low energy region ($<$10 keV). However, the detector response model built for ERs in XENON1T is, in principle, not applicable to $\gamma$-induced ERs that at sufficiently high energies may interact with the inner-shell electrons. When this happens, the vacancy in the inner shell results in either X-ray or Auger electrons emission, both of which further ionize xenon atoms. Consequently, $\gamma$-induced ERs can have multiple recoiling electrons instead of one as in $\beta$-induced ERs. The binding energy for L-shell electron in xenon is about 4.8-5.5 keV. According to the NIST database Wagner et al. (2003), the corresponding X-ray has mean free path of about 5 $\micro$m. The effect of the separation of electron clouds at this spatial scale on the recombination is not yet understood.

II.2 Detector Reconstruction Effects

Besides the intrinsic response of LXe, detector reconstruction effects on the S1 and S2 signals are modeled. More specifically, the spatial dependence of S1 and S2 signals, the single and double photoelectron (PE) emission of the PMT photocathode Faham et al. (2015); Paredes et al. (2018), the position reconstruction uncertainty, the reconstruction efficiency, bias, and signal fluctuations, and the acceptance of data selections in analysis are taken into account in the model.

Photons from an energy deposition and the subsequent recombination (Eq. 3) are detected by the PMTs as an S1 with an efficiency, which is the product of the light collection efficiency $\epsilon_{L}$, PMTs’ average quantum efficiency $\epsilon_{\mathrm{QE}}$, and PMTs’ average collection efficiency $\epsilon_{\mathrm{CE}}$. Electrons are drifted to the gas-liquid interface under the drift field, and are extracted under the stronger field, amplifying the electron signal (S2) by the gas gain $G$, which is the number of photoelectrons per electron that is extracted into gaseous xenon. Both $\epsilon_{L}$ and $G$ are spatially dependent, and related to the energy scale parameters $g^{\prime}_{1}$ (probability of one emitted photon to be detected as one PE) and $g^{\prime}_{2}$ (amplification factor for charge signal), respectively, by

[TABLE]

where $p_{\mathrm{dpe}}$ is the probability for the PMT photocathode to emit two photoelectrons when absorbing one photon Faham et al. (2015); Paredes et al. (2018), and $\epsilon_{\mathrm{ext}}$ is the extraction efficiency of the drifted electrons which is assumed to be constant in this study. Note that $g_{1}$ and $g_{2}$ in Aprile et al. (2019a) correspond to the averages of $g^{\prime}_{1}(x,y,z)$ and $g^{\prime}_{2}(x,y)$, respectively, in Eq. 10 over the active volume. The number of hits detected by PMTs, $N_{\mathrm{hit}}$, and photoelectrons generated from the PMT photocathode, $N_{\mathrm{pe}}$, can be described by a binomial distribution

[TABLE]

In addition to the ($x$, $y$) dependence caused by the varying charge amplification, S2 signals are a function of the $z$ position, because the electrons attach to electronegative impurities when drifting towards the gaseous phase. The number of electrons $N_{\mathrm{ext}}$ that survive the drifting and the extraction into the gas can be modeled as

[TABLE]

where $\tau_{e}$ and $\nu_{d}$ are the electron lifetime and electron drift velocity, respectively. The total proportional scintillation light detected, $N_{\mathrm{prop}}$, can be approximated as

[TABLE]

where $\Delta G$ is the spread of the gas gain. For simplicity, we consider $\Delta G/G$ a constant in the model.

The S1 and S2 signals are constructed from $N_{\mathrm{pe}}$ and $N_{\mathrm{prop}}$, respectively, amplified by the PMTs, digitized, and selected by clustering and classification software XENON Collaboration (2018). To account for biases and fluctuations in this process, the S1 and S2 are written as

[TABLE]

where $\delta_{\mathrm{s1}}$ ($\delta_{\mathrm{s2}}$) and $\Delta\delta_{\mathrm{s1}}$ ($\Delta\delta_{\mathrm{s2}}$) are the bias and spread, respectively, of the S1 (S2) reconstruction. Reconstruction biases and fluctuations are estimated using a waveform simulation including models of realistic scintillation pulse shape, charge amplification, electronic noise level, PMT single PE spectrum, PMT after-pulses, as well as secondary S2s induced by photonionizations on grids and impurities in the LXe volume XENON Collaboration (2018); Aprile et al. (2019a).

The S1 and S2 signals are corrected for their spatial dependence based on the reconstructed positions $\vec{x}_{r}$. The $z$ position of an event is reconstructed using the time difference between the S1 and the S2 signals, and has better resolution than the ($x$, $y$) position which is reconstructed through the S2 hit pattern on the PMTs in the top array. We assume the reconstruction fluctuations along $x$ and $y$ axes to be identical. The reconstructed position of $\vec{x}_{r}$ can be written as

[TABLE]

where $\sigma_{p}$ is the position reconstruction resolution, and depends on both S2 area and the ($x,y,z$) position of event. The corrected S1 (cS1) and S2 (cS2) are written as

[TABLE]

where $\tau_{e}^{\prime}$ is the mean electron lifetime measured. The uncertainty of the measured electron lifetime is used to constrain $\tau_{e}$ in the signal response model. The correction is based on the measured electron lifetime $\tau_{e}^{\prime}$. In the following analysis, as well as in XENON1T results Aprile et al. (2017b, 2018a), we use the corrected S2 collected by the bottom PMTs cS2b, which has a more homogeneous spatial dependence. The effect of field distortion is negligible and not implemented in the signal response model because the position correction to account for it is applied to data Aprile et al. (2018a).

Selection criteria were applied to data to ensure good quality of the sample and to optimize the signal-to-background ratio for the dark matter search. More details of the data selection can be found in Aprile et al. (2019a). The detection efficiency loss arises from the software reconstruction efficiency of S1s, and the S1-related and S2-related event selections. The efficiencies for these are considered as functions of $N_{\mathrm{hit}}$, S1 and S2, respectively, in the signal response model. In addition, a realistic selection of single scatters is implemented in the simulation. The rejection of multiple scatters is critical to the search for WIMP signals in XENON1T detector and is based on the areas of the largest and second largest S2s. Not all multiple scatters are rejected by this selection. This is mainly because of PMT after-pulses, photoionization of impurities in the detector, and the spatial resolution of the detector. Each energy deposition that is resolvable by the position reconstruction is taken into account in the simulation of the response model. The same single-scatter selection is applied to the simulated and actual data, in order to accurately address the acceptance of single scatters and the rejection power against multiple scatterings in the response model.

II.3 Fit to Calibration Data

The detector response model is constrained using the calibration data from 220Rn for ER and 241AmBe and a D-D generator for NR. Events with cS1 ranging from 0 to 100 PE are used to constrain the signal response model, covering the cS1 region of interest (3-70 PE) for WIMP searches Aprile et al. (2017b, 2018a). The detector response model obtained from the fit to calibration data is used to construct WIMP signal and background models, which are then input to the statistical inference of dark matter search data Aprile et al. (2017b, 2018a). The fit is performed simultaneously using all available XENON1T calibration data taken during the first (SR0 with drift field of 120 V/cm) and second (SR1 with drift field of 81 V/cm) science data taking periods. The fit uses the binned likelihood for distributions in log10(cS2b/cS1) versus cS1 using the data in the cylindrical fiducial volume defined in SR0 Aprile et al. (2017b). The likelihood is sampled using affine invariant Markov Chain Monte Carlo (MCMC) Goodman and Weare (2010). Important nuisance parameters in the detector response model are listed in Table 2. There are three parameters for scaling the S1 cut acceptance, S2 cut acceptance, and reconstruction efficiency. These parameters are constrained by the uncertainties estimated for the three efficiencies, which depend on the signal size. The $\delta_{\mathrm{s1}}$, $\delta_{\mathrm{s2}}$, $\Delta\delta_{\mathrm{s1}}$, and $\Delta\delta_{\mathrm{s2}}$ in Eq. 14 are signal-size dependent, and thus are not listed in Table 2. Their priors can be found in Aprile et al. (2019a). We choose to use the more conservative uncertainty between the lower and upper uncertainties given in Lenardo et al. (2015) for the NEST parameters that describe the response of LXe to NR, except for $\eta$, which parameterizes the Penning quenching of high-energy NRs. Correlations between the NEST parameters are not provided in Lenardo et al. (2015) and are, thus, not considered in the priors of this work, in order to be conservative and avoid potential over-constraint on fit parameters.

Figures  2 and  3 show the comparisons of the cS1 spectra and cS2b spectra, respectively, between the posterior of tested signal response model and data. For neutron calibrations the ER contamination are considered in the fit. For both ER and NR calibrations, a small fraction of events arise from accidental coincidence (AC) which will be illustrated in section III. The rates of each component are free in the fit, and are not listed in Table 2. The matching between the signal response model and the calibration data is good, with the goodness of fit (GoF) p-values (calculated using method in Gelman et al. (1996)) for the cS2b distributions comparison in different cS1 ranges shown in Table 1. Good agreement between model and data is critical for the WIMP signal and background models, especially for the (cS1, cS2b) region where the WIMP signal is expected. We also show the p-values for the match of the model with 220Rn in a reference (cS1, cS2b) region. The upper cS2b boundary of the reference region for matching 220Rn data is defined by the $97.7\%$ percentiles (+2$\sigma$) of NR events in 241AmBe calibration data. The GoF p-values for 241AmBe and D-D generator data are calculated excluding the lowest S2 region, corresponding to the $0.13\%$ percentile (-3$\sigma$) in S2. The GoF p-values at different cS1 for matching calibration data are all above the 5% threshold set for an acceptable fit.

III Background and Dark Matter Signal Models

Background and DM signal models are crucial in the statistical interpretation of the DM search results in XENON1T. There are four background components in XENON1T: ER, NR, surface, and accidental coincidence (AC). The ER and NR background models, as well as the WIMP signal model, are constructed based on the detector response model illustrated in section II (we call the detector response model “posterior” after fitting to the calibration data). The surface and AC background models are constructed using data-driven methods. The background and WIMP signal models are 3-D distributions in cS1, cS2b, and the spatial coordinate of the detector. The background and DM signal models also include the uncertainties in cS1, cS2b and spatial distribution, as well as in the absolute rate of the background and DM signal. In this section, the details of the background and DM signal models used in the statistical inference of XENON1T results Aprile et al. (2017b, 2018a) are given.

The final dark matter search is performed between 3$<$cS1$<$70 PE, and 50.1$<$cS2b$<$7940 PE. Below 3 PE, the S1 acceptance is very small due to the 3-fold PMT coincidence requirement for S1s. The upper $\mathrm{cS1}$ cut is chosen to contain most of spin-independent WIMP recoil spectra, shown in Fig. 8. In previous XENON analyses, fiducial volumes in radius R and $z$ were constructed to provide a low background for the analysis. With the inclusion of a model for the surface background, presented in Section III.3, the analysis could be extended to also consider the radius as an analysis variable. The analysis volume is defined by a maximal radius, $42.8~{}$cm, and a R-dependent $z$-cut, shown in Fig. 5 with a magenta line. The construction of this cut, which was made to include regions of the detector where the total background rate was approximately constant with $z$, is presented in Aprile et al. (2019a).

The magnitude of the radiogenic background, discussed in Section III.2, is attenuated moving towards the center of the detector. In order to optimize the discovery power of the analysis, a partition of the detector, with a clean, “core” volume was proposed. Optimizing for discovery significance yielded a central $0.65$ tonne volume, shown with a dashed green line in Fig. 5. The expected radiogenic neutron rate in this volume is $36\%$ of the average rate in the analysis volume.

III.1 Electronic Recoil Background Model

Although XENON1T achieved an excellent discrimination power between the ER background and NR signal, with an average ER leakage fraction below the NR median of about 0.3% Aprile et al. (2018a), the ER component is the dominant background for the DM search due to its high rate in comparison with the other background sources.

In the energy region of interest for WIMP search ($<$100 keVNR), the dominant component contributing to the ER background are $\beta$-decays of 214Pb. The 214Pb is a progeny of 222Rn, which is emanated from 238U daughters in the detector materials, and can convect and diffuse into the inner volume of the detector. Decays of 218Po and 214Bi-214Po, which are also 222Rn progenies, can be identified to estimate the rate of 214Pb. This selection is based on the unique energy and time profiles of 218Po and 214Bi-214Po decays, and gives activities of 71$\pm$5(stat.)$\pm$7(sys.) and 29$\pm$3(stat.)$\pm$3(sys.) events/ton/year/keV (tyu), respectively. The difference between the activities of 222Rn progenies is likely due to their plate-out onto the electrode and polytetrafluoroethylene (PTFE) reflector surfaces. As the 214Pb decay occurs between the 218Po and 214Bi decays, the rate of 214Pb is in range of 29 to 71 tyu which is consistent with the estimate of 56$\pm$6(stat.)$\pm$6(sys.) tyu from Aprile et al. (2016b), where 10 $\micro$Bq/kg of 222Rn were assumed.

The second largest component of the ER background is $\beta$-decays of 85Kr. The concentration of natural krypton natKr/Xe was reduced to 0.36$\pm$0.06 ppt by the end of SR0 through cryogenic distillation Aprile et al. (2017d). With regular rare-gas mass spectrometry measurements Lindemann and Simgen (2014) during SR1, we measured an average natural krypton concentration natKr/Xe of 0.66$\pm$0.11 ppt, resulting in an average decay rate for 85Kr of 7.7$\pm$1.3 tyu, using the conversion derived from data with high concentration of krypton at the beginning of XENON1T operation. These high-krypton concentration data also gave a 85Kr/natKr ratio of (1.7$\pm$0.3)$\times$10*-11* mol/mol. Taking the ER contributions from material radioactivity, solar neutrino scatterings, and $\beta$-decays of 136Xe of 8$\pm$1, 2.5$\pm$0.1, and 0.8$\pm$0.1 tyu, respectively, into account, the total ER background rate in the region of interest (ROI) for DM searches is estimated to be between 48$\pm$5 and 90$\pm$8 tyu. This is consistent with the prediction of 75$\pm$6 (sys.) tyu from Aprile et al. (2016b) and with the best-fit of 82${}_{-3}^{+5}$(sys)$\pm$3(stat) tyu of low energy ER background from Aprile et al. (2018a).

The energy distribution of the ER background in the ROI is assumed to be uniform due to the dominance of the flat spectrum from 214Pb beta decay. Uncertainties in the (cS1, cS2b) distribution for the ER background are dominated by the uncertainties in $\langle r\rangle$, and its fluctuations, $\Delta r$. The uncertainty in $\langle r\rangle$ is mainly from parameter $\gamma_{\mathrm{er}}$ in Eq. 8. The effects of varying $\gamma_{\mathrm{er}}$ and $\Delta r$ are shifting the mean and changing the spread, respectively, of the ER distribution in cS2b. Fig. 4 shows the variation of the ER distributions on log10(cS2b) in different cS1 ranges. The distributions are produced based on the detector response model with the rest of the nuisance parameters (shown in Table 2) marginalized to the point estimation (median posterior). Due to the computational complexity, the variation of background model in terms of (cS1, cS2b) distribution in the statistical interpretation is practically interpolated using the distributions that are computed at 2.3%, 6.7%, 15.9%, 30.9%, 50.0%, 69.1%, 84.1%, 93.3%, 97.7% percentiles of the posterior, which correspond to -2, -1.5, -1, -0.5, 0, 0.5, 1.0, 1.5, 2.0 $\sigma$, respectively, in standard deviations. Given that the ER background induced by radioactivities in detector materials is subdominant, the ER background is assumed to be spatially uniform inside the analysis volume.

III.2 Nuclear Recoil Background Model

The NR background, which has a similar (cS1, cS2b) distribution to the WIMP signal, contributes 1.43 events to the 1 tonne$\times$year exposure data Aprile et al. (2018a). Radiogenic neutrons, muon-induced neutrons and solar neutrinos contribute to this NR background.

Radiogenic neutrons are generated by ($\alpha$, n) reactions and spontaneous fissions of material radioactive impurities, and are the largest source of NR background. The neutron yields of materials are predicted using SOURCES-4A Wilson et al. (1999) based on the measured radioactivity of the detector materials Aprile et al. (2017a). The propagation of generated neutrons is simulated using the GEANT4 toolkit Agostinelli (2003).

Because of the large uncertainty ($\sim$50%) in the estimated neutron rate, multiple neutron scatter events in DM search data and calibration data are used to further constrain the rate uncertainty. For this purpose, multiple scatter events are unblinded prior to single scatters. Nine neutron multiple scatter events were identified in the DM search data (SR0 and SR1 combined) within the 1.3-tonne fiducial volume (FV). After a study comparing the single-to-multiple scatter ratio in data and simulations, the data-constrained neutron rate is estimated to be 55${}_{-20}^{+33}$ neutron scatters, including single and multiple scatters, per year in the 2-tonne active volume of XENON1T detector (the expectation is 37 neutron scatters per year from simulation). About 35% of the selected single scatters by radiogenic neutrons are misidentified multiple scatters. Fig. 5 shows the spatial and (cS1, cS2b) distributions of the identified neutron multiple scatters in DM search data together with the single scatter events surviving from the blinding cut Aprile et al. (2018a).

Neutrons emitted from bottom PMTs in the bottom array have a probability to scatter in the region between the TPC cathode and PMTs, referred to as the “below-cathode” region. The scintillation light from these scatters is detected, but electrons are lost because the electric field in the below-cathode region drifts them away from the active volume. Such events, with at least one scatter in the below-cathode region and a single scatter in the 1.3 tonne FV, are named neutron-X events, and have a lower cS2b to cS1 ratio than standard neutron scatters. The relative ratio and difference in (cS1, cS2b) distributions, as well as in spatial distributions, between normal neutron scatters and neutron-X events is taken into account in the simulation for building the NR background model.

Muon-induced neutrons are estimated to be subdominant to radiogenic neutrons in the active volume Aprile et al. (2016b), and the rate of muon-induced neutrons is further suppressed by applying a muon-veto selection. The final rate is approximately 2 orders of magnitude smaller than that of the radiogenic neutrons, and its contribution is neglected in the NR background modelling.

Compared with Aprile et al. (2017b), the NR background model induced by solar neutrino scatters has been updated with the latest results from COHERENT Akimov et al. (2017), and is 23% smaller. On top of the neutrino flux uncertainty of 14%, the updated NR background model also takes into account the uncertainty of about 21% from Akimov et al. (2017) and of about 22% from the signal response model at the energy region of interest for neutrino scatters. The resulting rate uncertainty for solar neutrino induced coherent elastic neutrino-nucleus scatters (CE$\nu$NS) is 34%.

III.3 Surface Background Model

Interactions in the bulk region of the detector can be modelled combining knowledge of LXe and detector-related responses, as discussed in the previous section for the ER and NR background models. For background events where the knowledge of the response is incomplete or missing, data-driven methods have to be developed to model the events. In this and the following section, two classes of such background components are analyzed for the XENON1T WIMP searches, events originating from the detector surface, and accidental coincidences of unrelated S1 and S2 signals.

Several experiments Agnese et al. (2015b); Amaudruz et al. (2018) have demonstrated that detector surfaces exposed to ambient air during construction are contaminated by a large amount of radon progeny, in particular 210Pb. With a 22 y half-life, 210Pb decays at a constant rate within the lifetime of the XENON1T experiment. For WIMP searches, ion recoils of 206Pb from 210Po $\alpha$-decays, $\beta$-decays and the resulting X-rays and Auger electrons of 210Pb are particularly important. Due to incomplete knowledge of LXe responses and detector physics in presence of PTFE, as well as complicated decay structures, a full model including the relevant physics processes has not yet been achieved in XENON1T. Instead, a data driven approach is adopted to predict the distribution of this background.

Background from radon progeny was modelled in cS1, cS2b, S2, R, and $z$ spaces, with a distribution $f_{\rm{Surf}}(\mathrm{cS1,S2,cS2}_{b},\mathrm{R},z)$. Radial position of surface events are reconstructed nearly symmetrically around the TPC boundary, with an uncertainty determined by the S2. Events mis-reconstructed outside the TPC are used to model the background distribution in S2, cS2b, cS1, and $z$, denoted as $f_{\mathrm{Surf-1}}(\mathrm{cS2}_{b},\mathrm{S2,cS1},z)$, with a kernel-density-estimation method Pedregosa et al. (2011). The distribution in cS2b and cS1 is shown in Fig. 6. Due to significant charge losses at the PTFE panels, the surface background overlaps significantly with the nuclear recoil region of interest (between nuclear recoil median and $-2\sigma$ quantile). In contrast, the R distribution provides excellent rejection power. To construct the distribution of surface background in R and S2 space, events originating at the PTFE surface are selected as control sample with an S1 size out of the region of interest, including the 210Po $\alpha$-decays. In each S2 slice, the radial distribution of the control-sample events is fitted, including an uncertainty estimated by different fitting functions. The 2-dimensional distribution of surface background, denoted as $f_{\mathrm{Surf-2}}(\mathrm{R,S2})$, is combined with $f_{\mathrm{Surf-1}}(\mathrm{cS2}_{b},\mathrm{S2,cS1},z)$ to form a complete model of the surface background as $f_{\mathrm{Surf}}(\mathrm{cS1,S2,cS2}_{b},\mathrm{R},z)$, including the uncertainty in the radial distribution.

The total rate is normalized to the number of events reconstructed outside the TPC boundary. In the WIMP search Aprile et al. (2018a), a radius cut of 42.8 cm is placed to reduce the surface background to $\sim 100$ events. In the likelihood fit described in section IV, the total surface background expectation value is conservatively treated as a free parameter.

III.4 Accidental Coincidence Background Model

The accidental coincidence (AC) of uncorrelated S1s and S2s (referred to as lone-S1 and lone-S2, respectively) is the fourth background component considered in the XENON1T model. Lone-S1s and -S2s originate from energy depositions in the non-active regions of the detector, where the scintillation or ionization signal are not detectable. For example, an energy deposition in the below-cathode region does not produce an ionization signal, and a deposited energy very close to the gate may have its scintillation signal blocked by the mesh. Similar to the surface background model described above, AC background is constructed through a data-driven approach by random pile-up of lone-S1 and lone-S2 samples from data.

The lone-S1 sample is obtained by searching for S1s in the time window before the larger primary S1 in each digitized event waveform. The estimated lone-S1 rate ranges from 0.7 to 1.1 Hz, depending on the requirement of noise rate in the search window. The difference in lone-S1 rate between science runs is negligible. The lone-S2 sample is obtained using events with no S1 found in the digitized waveform, or with a reconstructed $z$ position larger than the maximum drift time. The observed lone-S2 rate (with S2 threshold of 100 PE) is determined as 2.6 mHz and stays constant during both science runs. Selection criteria are applied in the lone-S1 and lone-S2 samples directly, while the event selections involving S2 (S1) are excluded. The AC event rate is calculated as $\mathrm{R_{AC}=R_{lS1}\cdot R_{lS2}}\cdot\Delta t$, where $\mathrm{R_{lS1}}$ and $\mathrm{R_{lS2}}$ are the rates of lone-S1s and lone-S2s, respectively. $\Delta t$ is the coincidence time window, which is fixed to the maximum drift time in the TPC (674 $\micro$s for SR0 and 727$~{}\micro$s for SR1). The AC background is nearly uniformly distributed in ($x$, $y$, $z$) position space. This yields a total AC background of 44.2 events in SR1, which is reduced to 7.0 events with an S2 threshold of 200 PE. Several selection criteria were developed to suppress AC background to a sub-dominant contribution and their rejection efficiency were estimated by simulation.

We simulate the distribution of AC events by sampling and randomly pairing the lone-S1 and -S2 samples. Interaction positions of the event are calculated by applying a field distortion correction to the sampled position. The correction of S1 and S2, depending on the event position, as well as event selections are all applied to the simulated sample. Examples are the drift-time dependent S2 width and S1 fraction in top PMT array cuts. The final AC background prediction is done with the simulation sample where all selection criteria are applied. The AC model was validated before unblinding, using both 220Rn calibration data and WIMP search data outside the ROI (for example the sample with S2 between 100 and 200 PE, or with S1s being identified as single electron S2s). The predicted AC rate after all selection criteria in the ROI for SR0 and SR1 combined is 0.47 - 0.74 events per tonne per year. As shown in Fig. 7, the AC distribution is concentrated in the low cS1 region, making it important in the search for light WIMPs.

III.5 WIMP Signal Model

The WIMP signal model depends on the dark matter mass, and assumes a uniform distribution of WIMP signals in the FV. The signal energy spectrum is computed in the same way as in Aprile et al. (2017b, 2018a). Distributions in ($\mathrm{cS1}$, $\mathrm{cS2_{b}}$), as well as the energy spectra for 10, 50, and 200 GeV/c2 WIMPs are shown in Fig. 8. The uncertainty in the WIMP signal model in the (cS1, cS2b) distribution, which comes from the uncertainty in the detector response model, is sub-dominant to the uncertainties of background models in the statistical inference. Therefore, we approximate the uncertainty of the WIMP signal model only in the form of its rate uncertainty. The approximated rate uncertainties as a function of WIMP mass are shown in Fig. 8.

IV Inference of XENON1T Data

In this section we describe the general techniques used in XENON1T for hypothesis testing and construction of confidence intervals. The focus is on the description of the specific likelihood used for the statistical interpretation of the XENON1T 1 tonne-year WIMPs search data Aprile et al. (2018a, 2019b).

IV.1 Hypothesis Testing and Confidence Intervals

The profile log-likelihood ratio is used as the test statistic for both confidence intervals and discovery assessment,

[TABLE]

where $\mathcal{L}$ stands for the XENON1T likelihood defined in Section IV.2, $\hat{\sigma}$ and $\hat{\boldsymbol{\theta}}$ are the signal and nuisance parameters that maximize the likelihood overall, and $\hat{\hat{\boldsymbol{\theta}}}$ are the nuisance parameters that maximize the likelihood with the condition that the signal strength is $\sigma$. To avoid unphysical regions, the best-fit signal $\hat{\sigma}$ is constrained to be non-negative.

A signal hypothesis $H_{\sigma}$ is tested against the data by computing the p-value of the observed test statistic $q(\sigma)$ given $H_{\sigma}$. In case of relatively large expected signal-like background, asymptotic formulae Cowan et al. (2011a); Wilks (1938) for the distribution of $q(\sigma)$ are convenient. However, given the low background of XENON1T, it was found that these approximations no longer hold. This led to an under-coverage of cross-section and resulted in a overestimated limit of $10\%$ on average in the SR0 results Aprile et al. (2017b). Therefore, the distribution of $q(\sigma)$ is computed using toy Monte Carlo (toy-MC) simulations of the background and signal models. Note that in computing these distributions the auxiliary measurements related to each nuisance parameter are also varied per toy-MC data sample.

A Feldman-Cousins construction in the profile likelihood Feldman and Cousins (1998); Patrignani et al. (2016) (also termed profile construction) is used to construct confidence intervals, using $q(\sigma)$. Sensitivities, as well as the coverage properties of the confidence band, detailed in Section IV.4, are explored with toy-MC simulations. Fig. 9 shows the distributions of lower and upper limits for a background-only simulation, illustrating the sensitivity computation. Both the science and calibration data-sets are drawn from their model distributions, as are the ancillary measurements according to their uncertainties. For the sensitivity computation, 1000 toy-MCs are run per mass point.

Before unblinding, the XENON1T collaboration resolved to use a higher threshold for reporting upper limits than the 1.28 $\sigma$ corresponding to a $90\%$ confidence level. Only the upper edge of the confidence interval is reported until the discovery significance reaches 3 $\sigma$. This leads to over-coverage at low WIMP cross section values, with an example shown in Fig. 11 affecting limits below the 1 $\sigma$ sensitivity band.

The FC construction may produce upper limits excluding signal strengths to which the experiment has a very small discovery power. In the case of a downwards fluctuation with respect to the background model, a $15\%$ power-constraint would be applied to set a lower threshold for the upper limit, as proposed in Cowan et al. (2011b).

The discovery significance is computed using the test-statistic evaluated at the null-hypothesis, $q_{0}\equiv q(\sigma=0)$. Denoting the distribution of $q_{0}$ under the null-hypothesis $H_{0}$ as $f_{H_{0}}(q_{0})$, the significance of a given observed test statistic $q^{\mathrm{obs}}_{0}$ can be expressed with a p-value,

[TABLE]

In practice, the p-value is estimated using toy-MC samples to account for non-asymptoticity due to the low signal-background overlap. For the FC construction, the test statistic distribution is estimated for 20 steps in the true signal for each mass. For each signal step, 2500 toy-MC simulations, including ancillary measurements are generated and fitted. A threshold curve, constructed by a smooth interpolation between the 90th percentile of the test statistic for each signal size is compared with the final log-likelihood ratio to construct the FC intervals.

The local discovery significance in Eq. 18 is computed for a single signal hypothesis. In the spin-independent analysis, the signal hypotheses are the WIMP masses considered between $6~{}\mathrm{GeV}/\mathrm{c}{}^{2}$ and $10000~{}\mathrm{GeV}/\mathrm{c}{}^{2}$. To compute a global discovery significance, the distribution of the most significant p-value for any hypothesis is estimated with 10000 toy-MC simulations, and compared with the result from data.

IV.2 The XENON1T Likelihood

The log-likelihood used in the spin-independent analysis is a sum of extended un-binned log-likelihoods for the two science runs. Additional terms are the extended un-binned likelihoods for ER calibration data, and terms expressing ancillary measurements of nuisance parameters $\theta_{m}$,

[TABLE]

where $\mathrm{SR}$ runs over data-taking periods, SR0 and SR1, $\sigma$ is the WIMP-nucleon cross-section, and $\mathcal{L}_{\mathrm{sci}}$, $\mathcal{L}_{\mathrm{cal}}$ are the likelihood terms for the WIMP search data and 220Rn calibration data, respectively. Ancillary measurements of nuisance parameters $\theta_{m}$ are included in $\log\mathcal{L}_{m}(\theta_{m})$. The un-binned science likelihood is defined in three dimensions: $\mathrm{cS1}$, $\mathrm{cS2_{b}}$ and R, the radius of the reconstructed event. Each background and signal distributions are defined and normalized in this three-dimensional space. The un-binned likelihoods take the form

[TABLE]

Here, the index $i$ runs over events in the relevant science data-set, and $c$ runs over each signal or background component, with expectation value $\mu_{c}(\sigma,\boldsymbol{\theta})$, which may be a nuisance parameter or a function of nuisance parameters. The probability density functions (PDFs), $f_{c}(\vec{x}_{i}|\boldsymbol{\theta})$, for each component are functions of the analysis coordinates $\vec{x}_{i}=(\mathrm{cS1}_{i},\mathrm{cS2}_{b}\,{}_{i},\mathrm{R}_{i})$ and are evaluated for each event in the likelihood. The models for the science data measurements, both the WIMP signal model, and the ER, neutron, CE$\nu$NS, surface and AC backgrounds are described in detail in the previous sections. In the case of the calibration likelihood, which utilizes a smaller volume, 0$<\mathrm{R}<$36.94 cm, the likelihood uses only the $\mathrm{cS1}$ and $\mathrm{cS2}$ dimensions, but the structure of $\mathcal{L}_{\mathrm{SR}\ \mathrm{cal}}$ is otherwise identical to Eq. 21.

The extension of detector modelling to include R was made to improve the sensitivity to a 50 $\mathrm{GeV}/\mathrm{c}{}^{2}$ WIMP by $10\%$, with respect to an optimized, smaller volume without radial modelling. The radial cut is set when the surface background starts dominating and no significant improvement in sensitivity can be obtained. In addition to the three analysis dimensions in which the distributions are modelled, events in the science data-sets are also classified as inside or outside a core mass of $0.65\ \mathrm{t}$ (shown in Fig. 5). Practically, this can be considered as a combination of two separate un-binned likelihoods in the shape of Eq. 21, where the relative expectations inside and outside the core mass are determined using the $\mathrm{R},z$ distribution of the different components.

IV.3 Nuisance Parameters

The signal and background models, consisting of expectation values and distributions in analysis space, depend on several nuisance parameters. Table 3 lists all the parameters of the combined likelihood, the data sets that mainly constrain them and their best-fit value. Uncertainties on the nuisance parameters are computed using the profiled likelihood in each nuisance parameter, but unlike the signal confidence intervals, the asymptotic construction is applied. Nuisance parameters are grouped in rate parameters, the expectation value for all background components, signal efficiency parameters and shape parameters that affect the model distributions. Expectation values for the 5 modelled backgrounds in the science data, as well as AC and ER rates in the 220Rn calibration likelihood, are all nuisance parameters in the likelihood. In the following, a separate nuisance parameter and term is applied for each science run, with the exception of the radiogenic rate.

The radiogenic rate, as well as the expected CE$\nu$NS rates for each science run, is constrained by ancillary measurements expressed as Gaussian likelihoods,

[TABLE]

These likelihoods are defined by the PDF for the ancillary measurement of the component $\hat{\mu_{c}}$, given a true expectation value $\mu_{c}$ and measurement uncertainty $\sigma_{c}$. The signal expectation is

[TABLE]

where $\mu_{\mathrm{ref}}$ is the reference expectation for a WIMP of mass $M$ and cross-section $\sigma_{\mathrm{ref}}$, given by the signal acceptance, and $\epsilon$ is a multiplicative factor expressing the uncertainty on the signal expectation for a fixed cross-section. This expectation uncertainty is constrained by the NR model posterior taking all model variations into account,

[TABLE]

where the uncertainty $\sigma_{\epsilon}(M)$ depends on the WIMP mass, ranging from $0.15$ at $6~{}\mathrm{GeV}/\mathrm{c}{}^{2}$ to $0.03$ at $200~{}\mathrm{GeV}/\mathrm{c}{}^{2}$, as shown in the inset of Fig. 8. The nominal value of $1$ reflects that the best-fit expectation is expressed by the reference expectation.

The AC rate $\mu_{\mathrm{AC}}$ is constrained between two extreme estimates of its rate and is assigned a uniform PDF between the lower and upper reference, written as $0.6\cdot\hat{\mu_{\mathrm{AC}}}$ and $\hat{\mu_{\mathrm{AC}}}$ for convenience,

[TABLE]

The ER and surface background rates are not assigned auxiliary measurements, as their high statistics in the science data sample constrain them. In the case of the surface background, the region of highest signal overlap (at low $\mathrm{cS1}$ and R) was blinded. This motivated the conservative procedure of not placing an auxiliary constraint on the surface shape.

In addition to uncertainties on the rate, the ER and surface background PDFs are also assigned shape uncertainties. The surface background sideband fit, described in Section III.3, is also used to construct an uncertainty for the radial slope of the background. To avoid over-constraining this distribution, which could lead to spurious excesses or too-tight confidence intervals, no constraint is placed on the radial slope from the sideband measurement. The science data fit finally constrains the radial slope with approximately $2.5$ times smaller uncertainty than the sideband fit.

The ER model is described in detail in Section III.1. The nuisance parameters $\gamma_{\mathrm{er}}$ and $\Delta r$ described there are propagated to the likelihood, with the former shared between science runs. The value of these nuisance parameters and their associated errors are determined in the combined fit by including the ER calibration likelihood term for 220Rn data. This ER model includes the nominal ER model, as well as variations due to changing the photon yield or recombination fluctuation parameter. It is slightly different from the one illustrated in Section II, as the uncertainties in the model from nuisance parameters other than $\gamma_{\mathrm{er}}$ and $\Delta r$ are not included in the nominal ER model.

A mismodelling term, or “safeguard”, proposed in Priel et al. (2017), is a shape-uncertainty added to the ER model, consisting of a signal-like component added to or subtracted from the nominal ER model. This ensures that regardless of other nuisance parameters, the ER model will have the freedom to fit the calibration data in the signal-like region. A spurious signal-like over- or under-estimation in the ER background will have the greatest impact on the inference, giving spuriously constraining limits or spuriously significant excesses, respectively. The safeguard mainly affects the ER model tail which overlaps with the signal region, as shown in Fig. 10. The ER model, including the mismodelling term is constrained by the calibration likelihood terms included in the total likelihood, this allows the safeguard component to be constrained by the much higher (approximately 10 times for SR0 and SR1) statistics of the 220Rn calibration data compared to the science data. The ER background PDF including the safeguard term, PDF $f^{\prime}_{\mathrm{ER}}(\vec{x}|\boldsymbol{\theta}_{\mathrm{ER}})$ can be written as

[TABLE]

where $a_{\mathrm{mm}}$ is the safeguard nuisance parameter and $f_{\mathrm{sig}}$ the signal PDF. If the safeguard causes the PDF to be negative in a region, it is truncated to 0. The pre-factor $b$, a function of $a_{\mathrm{mm}}$ and nuisance parameters that affect the ER distribution $\boldsymbol{\theta}_{\mathrm{ER}}$, ensures that the PDF is normalized in the analysis space.

The compatibility between the best-fit to data and the safeguard-equal-zero hypothesis is performed using the profiled log-likelihood ratio to compute the 1-sigma error on the best-fit safeguard. The calibration data indicates that the signal-like tail of the ER is less pronounced than in the nominal model, reflected in a negative safeguard fit. A zero safeguard is $2.9\,\sigma$ from the SR1 +SR0 combined fit using this method.

IV.4 Coverage

The fraction of repeated experiments where the confidence interval contains the true parameter is called the coverage. Perfect coverage is equal to the confidence interval, $0.9$ in the case of XENON1T. While the nominal FC construction, introduced in section IV.1, provides coverage by construction, the coverage of the profile construction must be investigated for the likelihood in question. As we decided to report only the upper edge of the confidence interval for discovery significances $<3\,\sigma$, there will be over-coverage at very low signal sizes. This has a similar effect to the power constraint. Fig. 11 shows the coverage for a $50$ $\mathrm{GeV}/\mathrm{c}{}^{2}$ WIMP, both for the profile construction (orange points), and the analysis including the $3\,\sigma$ threshold (blue points). The green band shows the $-1\,\sigma$ to $1\,\sigma$ sensitivity band. The result is consistent with perfect coverage (black line), with overcoverage for the $3\,\sigma$ threshold only under the $-1\,\sigma$ edge of the sensitivity band. The effect on coverage of mismeasuring the nuisance parameters was also studied. In the case of the mismodelling term, a shift in the parameter by three times the observed value was required for a one percentage point shift in the coverage.

V Summary

In this manuscript we have reported details of the detector response model, the background and WIMP signal models, and the statistical inference of XENON1T analysis chain. These have been used in the interpretation of results in the search for spin-independent elastic WIMP-nucleon interactions in XENON1T Aprile et al. (2017b, 2018a), and is also used in searches for alternative dark matter candidates or interactions using XENON1T data Aprile et al. (2019c, b). The background model, using a simulation-based detector response model, has been improved with respect to previous analyses by reduction of the systematic bias and a better treatment of parameter uncertainties and correlations. The statistical inference was developed to use a combined, unbinned likelihood that is adaptable to multiple analysis spaces and additional data sets. A signal-like mismodeling term is introduced for the first time in an analysis as a background model shape uncertainty. In addition, the inference now employs toy Monte Carlo simulations extensively to construct and validate the confidence bands.

Acknowledgement

We gratefully acknowledge support from the National Science Foundation, Swiss National Science Foundation, German Ministry for Education and Research, Max Planck Gesellschaft, Deutsche Forschungsgemeinschaft, Netherlands Organisation for Scientific Research (NWO), Netherlands eScience Center (NLeSC) with the support of the SURF Cooperative, Weizmann Institute of Science, Israeli Centers Of Research Excellence (I-CORE), Pazy-Vatat, Initial Training Network Invisibles (Marie Curie Actions, PITNGA-2011-289442), Fundacao para a Ciencia e a Tecnologia, Region des Pays de la Loire, Knut and Alice Wallenberg Foundation, Kavli Foundation, and Istituto Nazionale di Fisica Nucleare. Data processing is performed using infrastructures from the Open Science Grid and European Grid Initiative. We are grateful to Laboratori Nazionali del Gran Sasso for hosting and supporting the XENON project.