XENONnT WIMP Search: Signal & Background Modeling and Statistical Inference

Abstract

The XENONnT experiment searches for weakly-interacting massive particle (WIMP) dark matter scattering off a xenon nucleus. In particular, XENONnT uses a dual-phase time projection chamber with a 5.9-tonne liquid xenon target, detecting both scintillation and ionization signals to reconstruct the energy, position, and type of recoil. A blind search for nuclear recoil WIMPs with an exposure of 1.1 tonne-years yielded no signal excess over background expectations, from which competitive exclusion limits were derived on WIMP-nucleon elastic scatter cross sections, for WIMP masses ranging from 6 GeV/$c^{2}$ up to the TeV/$c^{2}$ scale. This work details the modeling and statistical methods employed in this search. By means of calibration data, we model the detector response, which is then used to derive background and signal models. The construction and validation of these models is discussed, alongside additional purely data-driven backgrounds. We also describe the statistical inference framework, including the definition of the likelihood function and the construction of confidence intervals.

I Introduction

Astrophysical and cosmological observations at various scales provide clear evidence for the existence of dark matter (DM) as a fundamental building block of our universe [1, 2]. Numerous extensions of the standard model of particle physics predict additional fundamental particles as potential candidates for DM [3, 4]. Measuring the direct interaction of particles from the DM halo of the Milky Way with a suitable detector could provide conclusive evidence for the particle DM hypothesis. The XENONnT experiment aims to detect a signal from weakly interacting massive particles (WIMPs) elastically scattering off xenon nuclei with a detector operated underground at the INFN Laboratori Nazionali del Gran Sasso (LNGS) in Italy. The collaboration published the first WIMP search results from a data-taking period between July $6^{\mathrm{th}}$ and November $10^{\mathrm{th}}$ 2021 with a total live time of $95.1$ days and a fiducial mass of $(4.18\pm 0.13)\;\mathrm{t}$, referred to as science run 0 (SR0). We performed a blind analysis and observed no significant excess above background expectations [5]. This led to a minimum upper limit on the spin-independent WIMP-nucleon cross section of $2.58\times 10^{-47}\;\mathrm{cm}^{2}$ for a WIMP mass of $28\text{\,}\mathrm{GeV}\text{/}{\mathrm{\text{$c$}}}^{2}$ at $90\%$ confidence level.

The detector used for the XENONnT experiment is a dual-phase xenon time projection chamber (TPC). A particle interaction in the active $5.9\text{\,}\mathrm{t}$ liquid xenon (LXe) target results in prompt vacuum ultraviolet scintillation light as well as free ionization electrons. The instantaneous light signal, called the S1 signal, is detected by two arrays of photomultiplier tubes (PMTs) at the top and bottom of the approximately cylindrical TPC. An electric field ($23\text{\,}\mathrm{V}\mathrm{/}\mathrm{c}\mathrm{m}$ average electric drift field during SR0) is applied to the target volume to drift the ionization electrons toward the liquid xenon surface. Here, they are extracted and accelerated into a xenon gas volume above the liquid via a stronger electric field ($2.9\text{\,}\mathrm{k}\mathrm{V}\mathrm{/}\mathrm{c}\mathrm{m}$ electric extraction field during SR0). The kinetic energy of the accelerated electrons in the gas phase is sufficient for the emission of electroluminescent light, which is proportional to the number of extracted electrons. This second light signal is referred to as the S2 signal. The distribution of the S2 signal detected by the PMTs in the top array is used to infer the position of the initial interaction in the horizontal plane $\mathrm{(\mathrm{X},\mathrm{Y})}$ parallel to the liquid surface, which defines the radial coordinate $\mathrm{\mathrm{R}=\sqrt{\mathrm{X}^{2}+\mathrm{Y}^{2}}}$. The vertical component $\mathrm{Z}$ follows from the drift time of electrons, which is determined as the time difference between the S1 and S2 signals. This three-dimensional position reconstruction of the interaction vertex enables the discrimination of multi-site events as well as the selection of an inner fiducial volume (FV) within the TPC. This volume benefits from a particularly low background level thanks to the excellent self-shielding properties of LXe, given its high density. Moreover, the relative magnitudes of the S1 and S2 signals can be used to determine if a particle interacted with the xenon nucleus (expected from WIMPs) or xenon shell electrons (typical of the dominant sources of background), i.e. nuclear recoil (NR) or electronic recoil (ER) events. The electric drift and extraction fields are established by three parallel-wire electrodes (cathode, gate, and anode from bottom to top). To reduce wire sagging, two and four horizontal perpendicular wires support the gate and anode electrodes, respectively.

The TPC is nested in an active neutron veto, which in turn is housed in an active muon veto. Both veto systems are water Cherenkov detectors. More details on the TPC, the veto systems, and other subsystems of the detector are provided in [6, 7].

The data analysis chain of the WIMP DM search in XENONnT is subdivided into two major parts, like in XENON1T [8, 9]. The first part is discussed in [10] and comprises signal and event reconstruction, corrections, event selection, and energy scale calibration. Here we present the second part of the analysis chain: The detector response to the different interaction types is modeled based on calibration data and discussed in Sec. II. The derived best-fit models are important inputs for the definition of signal and background models, detailed in Sec. III. This section also covers data-driven background models that do not rely on the detector response. Finally, the statistical methods used to compute constraints on DM given the XENONnT data are discussed in Sec. IV.

II Modeling the detector response

Modeling the detector response to energy depositions allows for a powerful discrimination between NR signal and ER background in S1–S2 space. NR events are produced by particles scattering elastically off xenon nuclei, while in ER events, particles scatter off xenon shell electrons. The different energy loss processes involved after the two types of recoil cause ER and NR events to form separate populations in S1–S2 space, which is exploited in the WIMP DM search.

In this section, we describe the TPC response model to energy depositions in LXe via ER or NR interactions, which is obtained from fits to calibration data. This detector response model is separated into two parts. The first one is the empirically parametrized LXe emission model, which describes the production processes of the detectable quanta, i.e. the conversion of deposited energy into scintillation photons and ionization electrons. The second part is the detector reconstruction model, which covers the conversion from the produced photons and electrons into the observed and spacetime-dependence corrected S1 and S2 signals (cS1 and cS2). Due to its complexity, it is unfeasible in the fit to simulate processes of photon and electron propagation on a per-quantum basis. Instead, toy-Monte-Carlo (toy-MC) simulations of the detector reconstruction model are used which sample from effective maps provided either by data-driven analyses described in [10] or by simulation-driven analyses using the XENONnT Monte Carlo (MC) framework [11, 12]. Then all the model parameters in the simulations are fit to both ER and NR calibration data. The two parts of the detector response model together with the toy-MC simulation workflow are described in Secs. II.1 and II.2, and the fit to ER and NR calibration data is detailed in Sec. II.3.

II.1 Liquid Xenon Emission Model

The production of quanta from energy depositions in a xenon target is complex and lacks a comprehensive description derived from first principles. Thus, a semi-empirical parametrization of the emission model is commonly used. The parametrization used in XENONnT SR0 is similar to XENON1T [9], which is based on the Noble Element Simulation Technique (NEST) model described in [13, 14]. The simulation workflow of the emission model is described in detail in Appendices A and B, and is summarized in the following.

In an ER event, the recoil energy is transmitted into the excitation and ionization of xenon atoms. Recoiling xenon nuclei from an NR event, on the other hand, lose a significant amount (roughly 80 %) of their kinetic energy to atomic motion in collisions with other xenon atoms [13]. This thermalization process is undetectable in LXe TPCs, resulting in an effective signal loss. For both recoil types, the total number of detectable quanta $N_{\mathrm{q}}$ equals the sum of the number of produced excitons $N_{\mathrm{ex}}$ and ions $N_{\mathrm{i}}$. A fraction of ions and electrons recombine, resulting in the production of additional excitons. An exciton can combine with a ground-state xenon atom forming an excimer, which de-excites, producing a scintillation photon. In an event, $N_{\gamma}$ photons are generated, part of which are then detected by PMTs as the $\mathrm{S1}$ signal. The $N_{\mathrm{e}}$ electrons that do not recombine are drifted to the liquid-gas interface and eventually form an $\mathrm{S2}$ signal.

In the region of interest (ROI) for the SR0 WIMP search, $\mathrm{cS1}$ in [0, 100] PE and $\mathrm{cS2}$ in [$10^{2.1}$, $10^{4.1}$] PE, NR interactions have a larger recombination probability than ER interactions given the same total energy deposition. This leads to smaller ratios of S2 to S1 areas, which is the principle behind ER-NR discrimination in LXe TPCs. We parametrize the mean recombination fraction for ER interactions using the Thomas-Imel box model [15] with an additional Fermi-Dirac term, as in [9]. For NR interactions, the emission model follows the parametrization used in the NEST v1 model detailed in [14]. While there are newer versions of fits available in the NEST framework that use a different parametrization for NR interactions (which we refer to as NEST v2 in the following), we choose to stay with the previous version (referred to as NEST v1). The model is simpler, fits our data well, and the best-fit values and uncertainties of parameters of NEST v2 were not published at the time of the analysis.

II.2 Detector Reconstruction Model

The detector reconstruction model covers all aspects of the signal reconstruction, starting from the produced scintillation photons and ionization electrons up to the measured S1 and S2 signal sizes. The XENONnT detector reconstruction model is almost identical to the one presented in [9] and is briefly outlined here. The detailed simulation workflow of the detector reconstruction model is described in Appendix C.

For S1 signals, the spatial dependence of geometric effects during photon propagation, the photon detection efficiency of the PMTs, and the effect of double photoelectron emission (DPE) from the PMT photocathode [16, 17] are all modeled in the simulations. For S2 signals, we model the loss of electrons along their drift path due to attachment to electronegative impurities as well as due to electric field effects close to the TPC wall, the efficiency of electron extraction from the liquid into the gas phase, and the single-electron gain of extracted electrons to detected photoelectrons (PE) via the electroluminescence process in gas. In addition to these physical processes, we also account for the influence of software reconstruction effects on the signals. For both signal types, the software reconstruction process can introduce biases and fluctuations in the recorded signal size. We also account for the S1 signal reconstruction efficiency, which vanishes at about 3 PE, as we require that at least 3 PMTs detect a photon within $\pm$50 ns around the maximal amplitude. Finally, we model the uncertainty of the event position reconstruction and the acceptances of data quality selections [10].

The modeling of ER interactions relies on calibration data from sources dissolved in xenon, in particular ²²⁰Rn and ³⁷Ar. The ²²⁰Rn progeny ²¹²Pb undergoes a $\beta$-decay with a Q-value at about 0.6 MeV [18], giving an approximately flat ER energy spectrum in the WIMP ROI, i.e. below about $10$ keV. This dataset is used for the fit of the ER model. The ³⁷Ar source provides a mono-energetic ER peak when its K-shell electron capture process leads to a fast cascade of X-rays and Auger-Meitner electrons, with a total energy of about 2.8 keV [19]. Both ER sources provide single-scatter (SS) events, producing one S1 and one S2 signal for the total deposited energy, and the above-described detector reconstruction processes are simulated accordingly.

The NR model is calibrated using an external ²⁴¹AmBe neutron source (referred to as AmBe source in the following). It emits a fast neutron via the ⁹Be$(\alpha,n)^{12}$C capture reaction. Neutrons often scatter multiple times in the LXe target in so-called multi-scatter (MS) events, creating distinct energy depositions at each interaction site, thus the topologies of neutron events must be considered in the modeling. Since the neutrons travel at about 10 % of the speed of light and the S1 width is $\sim$$\mathcal{O}(\mathrm{100\,ns})$, the S1 signals from the separate energy depositions in a MS event get reconstructed into one merged S1 signal. The S2 signals might be distinguished based on the separation in spacetime between the energy depositions of the individual scatters, but can get (partially) merged if the interaction sites are close together. Due to the low drift field of $23\,\mathrm{V/cm}$ in XENONnT SR0, the $\mathrm{Z}$ separation resolution between two S2 signals is reduced with respect to the design because of the increased drift time $\sim$$\mathcal{O}(\mathrm{ms})$ and consequently larger S2 spread $\sim$$\mathcal{O}(\mathrm{\upmu s})$ due to the diffusion of electrons. This effect is accounted for in the fit of the NR model to the neutron calibration data: we first produce the photons and electrons for each energy deposition in a MS event separately, based on the emission model. While the primary scintillation photons are summed to produce the merged S1 signal, we only sum S2 signals that will become part of the largest S2 signal after going through the remaining detector reconstruction processes. This information is provided by PMT waveform simulations of neutron scatters from the AmBe source, using the Monte Carlo (MC) framework of XENONnT [11, 12] (see also Sec. III.3 for details on neutron simulations). The results of these waveform simulations are used as inputs to the toy-MC simulations and fitting framework. Data selections that remove MS events (both resolved and unresolved) are directly applied to these inputs, because their acceptances can only be determined from full simulations, and are correspondingly dropped from the selection acceptance curves for the NR fit.

Figure 1: AmBe neutron calibration events with (red) and without (gray) a coincident signal in the neutron veto (NV). Selecting coincident events ensures a clean nuclear recoil sample for detector response modeling. The accidental coincidence population (cS1 below 5 PE) and misidentified single-electron S1s (cS1 between 10 and 30 PE) are visible in the gray population. After applying additional data quality selections, approximately 2000 events (Fig. 2, right) are used for band fitting.

II.3 Fit to Calibration Data

Figure 2: Comparison between calibration data and the best-fit ER (left) and NR (right) models. The equiprobable binning for the 2D binned Poisson likelihood $\chi^{2}$ goodness-of-fit tests is shown. The color scale indicates the deviation of the number of data points (overlaid as black dots) in each bin from the best-fit model expectation $\mu_{\mathrm{bin}}$ in units of sigma, which is 29.3 for the ER model and 36.9 for the NR model.

The parameters in the full detector response model are fitted to ²²⁰Rn, ³⁷Ar, and AmBe calibration data. The ROI selection is the same as the SR0 WIMP search region defined in Sec. II.1. All data quality selections are applied to the calibration data within the ROI. For the AmBe neutron data, an additional selection is applied using coincident gammas detected in the neutron veto. In ⁹Be$(\alpha,n)^{12}$C capture reactions in the AmBe source, ¹²C is left in the first excited state after neutron emission in about 50 % of the cases, estimated with neutron veto data [20]. The ¹²C excited state promptly decays via the emission of a 4.44 MeV gamma ray, coincident with the neutron emission. The requirement of the time coincidence within $250\,\mathrm{\upmu s}$ between the 4.44 MeV gamma ray detected in the neutron veto and the event in the TPC, leads to a clean calibration dataset, with only about one expected event from accidental coincidences between the two detectors in the dataset ($\sim$0.05 % of the resulting sample). The events recorded in the ROI that pass (fail) the neutron veto coincidence selection are shown in red (gray) in Fig. 1.

For the fit to the ER data, an additional background component is added to account for accidental coincidence (AC, see subsection III.5 for details) events formed by a wrong pairing of S1 and S2 signals in the calibration dataset.

After all selections, the ²²⁰Rn and AmBe datasets consist of $\sim$2000 events each. We use a downsampled dataset of 10 000 ³⁷Ar events to perform a combined ²²⁰Rn-³⁷Ar ER fit. With the downsampling we avoid overconstraining the fit to the narrow, very low energy range of ³⁷Ar.

Figure 3: Photon (top) and charge (bottom) yields as functions of deposited energy for ER events (left) and NR events (right). Dark red lines and shaded bands show the XENONnT liquid xenon emission models with the best-fit parameters and their uncertainties, respectively. The fit results from XENON1T (blue shaded bands) [9, 21], and the models from NEST v1 (gray solid line) and NEST v2 (black solid line), both evaluated at 23 V/cm, are shown as well. Several measurements from literature performed at various drift fields as labeled in the legends (ER yields from [9, 19, 22, 23, 24, 25], NR yields from Aprile 2005 [26], Aprile 2006 [27], Aprile 2009 [28], XENON100 [29], Plante 2011 [30], Sorensen 2009 [31], Manzur 2010 [32], LUX 2016 [33], LUX 2022 [34]) are also shown. The gray shaded regions mark the energy ranges in which the total detection + selection efficiency in the WIMP analysis is above 10 %.

For the fit, an affine-invariant Markov chain Monte Carlo (MCMC) algorithm [35, 36] is used to sample from the high-dimensional posterior distribution of the parameters. In each step of the sampling, a GPU-accelerated toy-MC simulation is performed for every set of the ER and NR model parameters, following the steps described in the previous two sections, producing a model of the detector response in the space of corrected signals, i.e. $\mathrm{cS1}$–$\mathrm{cS2}$. The model parameters are described in detail in the appendix and are listed in Tab. 2. These toy-MC models are then compared to calibration data by calculating the binned Poisson likelihood in $\mathrm{cS1}$–$\mathrm{cS2}$ space. The likelihoods are multiplied by the prior distributions of the parameters to yield the posteriors. After sufficient iterations of the MCMC, the samples in the chain then converge to the posterior distribution of the parameters.

In the ER emission model, the parametrization is the same as in XENON1T [9] with very wide flat priors, referred to as free priors. For the NR emission model, the parameter priors are taken from  [14]. A Gaussian-like distribution with different widths on either side is used as the prior of parameters with asymmetric uncertainties. However, due to the low drift field of 23 V/cm, where no literature data on NR yields exist, the validity of the field dependence in the model is unverified. Therefore, the field dependence in the emission model from [14] was modified. Two scaling parameters were freed ($\alpha$ and $\gamma$ in Eq. 24 and Eq. 27 in the appendix), and two field dependence parameters ($\zeta$ and $\delta$) were fixed to the reported best-fit values. The widths of all other parameter priors were doubled in order to allow more freedom in the fit. Because the dependence of NR signals on the drift field has been shown to be small (see also literature measurements shown in Fig. 3), we consider this a reasonable choice.

The fit results of all emission model parameters are summarized in Tab. 2 in the appendix. No significant tension between the posteriors and the priors of any parameter is found. Suitable goodness-of-fit (GOF) tests are performed to assess whether the best-fit models adequately describe the ²²⁰Rn and AmBe calibration data. Specifically, we employ binned Poisson likelihood $\chi^{2}$ tests in the $\mathrm{cS1}$–$\mathrm{cS2}$ space. We adopt an equiprobable binning scheme using the GOFevaluation package [37], ensuring that the number of expected events under the best-fit model is the same in each bin, $\mu_{\mathrm{bin}}$. In Fig. 2 the results of the 2D $\mathrm{cS1}$–$\mathrm{cS2}$ Poisson likelihood $\chi^{2}$ test are shown for both the ER and the NR fit, overlaid with the ²²⁰Rn and the AmBe calibration data, respectively. The resulting p-values of the tests indicate no significant discrepancy between the best-fit models and the calibration data.

Figure 3 shows the photon and charge yields as a function of deposited energy for the ER and NR emission models using our best-fit parameters (dark red). For comparison, NEST models calculated at a drift field of 23 V/cm and several published measurements at different drift fields are shown. For the NR yields, we also compare our results to the NEST v1 model at a drift field of 23 V/cm (red dashed lines in Fig. 3) which uses the same parametrization as the model in this work (see Sec. II.1), showing very good agreement.

Figure 4: Distribution of each background component (colored) in $\mathrm{cS1}$–$\mathrm{cS2}$ space. The NR background includes both radiogenic neutron background and NR events by neutrinos. The probability density functions (PDFs) of backgrounds are shown as 1$\sigma$ (2$\sigma$) dark (light) regions, containing 68% (95%) expected events in the ROI. For reference, the contours are also shown for a spin-independent 200 GeV$/c^{2}$ WIMP as a solid (dashed) dark gray line. The light gray dash-dotted lines are contours of constant NR-equivalent energy (and ER-equivalent energy in the top left panel) in units of $\mathrm{keV}$.

III Signal and background modeling

The dominant sources of background in the WIMP search are ER events from intrinsic $\beta$-decays from materials and neutrinos (Sec. III.2), NR events from radiogenic neutrons from detector materials (Sec. III.3) and coherent neutrino scattering (Sec. III.4), AC events (Sec. III.5), and surface background (Sec. III.6). Besides their total expected rates, the distribution of these background events in the analysis space cS1–cS2–$\mathrm{R}$ is derived. In Fig. 4, the distributions of background models in the $\mathrm{cS1}$–$\mathrm{cS2}$ space are shown, compared to the signal model for elastic spin-independent scattering of a 200 GeV$/c^{2}$ mass WIMP.

III.1 WIMP Signal Model

For non-relativistic, spin-independent elastic WIMP-nucleus coherent scattering, the event rate $R_{\mathrm{WIMP}}$ scales with the square of the atomic mass number $A^{2}$ of the target and can be written as [38]

$$\frac{\mathrm{d}R_{\mathrm{WIMP}}}{\mathrm{d}E}=\frac{\rho_{\chi}}{2m_{\chi}\mu_{\chi N}^{2}}\left\langle\frac{1}{v}\right\rangle\sigma A^{2}F^{2}(q),$$ (1)

where $E$ is the recoil energy, $\rho_{\chi}$ is the local DM density of 0.3 GeV/($c^{2}\,$cm³), $m_{\chi}$ is the WIMP mass, $\mu_{\chi N}$ is the WIMP-nucleon reduced mass, $v$ is the WIMP velocity in the lab-frame, $\sigma$ is the WIMP-nucleon cross section, and $F(q)$ is the nuclear form factor as a function of the momentum transfer $q$ to the xenon nucleus [39]. The DM velocity distribution is averaged using the parameter values of the standard halo model with values from [40, 41, 42, 43, 44, 45], as recommended in [46]. Fig. 5 shows the cS1–cS2 distribution of WIMP-nucleus scattering for different WIMP masses. For increasing masses, the $68\%$ contours extend to higher cS1 and cS2 values, up to about 200 GeV/$c^{2}$, where the shape does not change significantly anymore due to kinematic effects. For this reason, confidence intervals in $\sigma$ (see Sec. IV) for $m_{\chi}\gtrsim 200\,\mathrm{GeV}/c^{2}$ approximately scale proportional to the inverse of the assumed WIMP mass $m_{\chi}$.

For the spin-dependent WIMP-nucleus interaction, assuming natural abundances of xenon isotopes, only ¹²⁹Xe and ¹³¹Xe can contribute since they are the only two stable isotopes of xenon with non-zero nuclear spin $J$. The cross section is usually written as

$$\frac{\mathrm{d}\sigma_{\mathrm{SD}}}{\mathrm{d}q^{2}}=\frac{\sigma}{3\mu_{\chi N}^{2}v^{2}}\frac{\pi}{2J+1}S_{N}(q).$$ (2)

The axial-vector structure factor of xenon $S_{N}$ is taken from [47]. In the XENONnT SR0 analysis, the mean of the structure factor is used and the uncertainty is neglected since it is much smaller than the uncertainty from the NR fit.

Uncertainties from the posterior distribution of the NR model parameters and other efficiencies can be propagated to the NR rate uncertainty. For a 6 GeV$/c^{2}$ (50 GeV$/c^{2}$) WIMP, the relative rate uncertainty is $\sim$30% (10%). The NR model shape can also be affected by the posterior distribution. However, because of the low statistics of NR events in SR0 (for both background and signal expectation), the impact of the shape uncertainty of the NR model is negligible compared to the uncertainty of its absolute rate. Thus, the NR model shape uncertainty is not propagated to the final likelihood.

Figure 5: Expected distribution in cS1–cS2 space from spin-independent WIMP nuclear recoil for WIMP masses from 6 GeV/$c^{2}$ to 1000 GeV/$c^{2}$. The contours contain 68% of expected events in the ROI. The blinded region is shown for reference as a gray band and dash-dotted lines are constant NR-equivalent energy contours.

III.2 Electronic Recoil Background Model

Figure 6: Shape uncertainty of the ER model constrained by the ²²⁰Rn calibration data. The blue (gray) solid and dashed lines represent the contours of ER background (200 GeV$/c^{2}$ WIMP) containing 68% and 95% of expected events in the ROI. The transparent blue regions show the shape uncertainty by varying the two shape nuisance parameters by $\pm 1\sigma$. Dash-dotted lines are contours of constant ER-equivalent energy in units of $\mathrm{keV}_{\mathrm{ER}}$.

ER events are one of the dominant background sources in the XENONnT SR0 WIMP search, primarily originating from $\beta$-decays of intrinsic radioactive contaminants such as ²¹⁴Pb (a product of the ²²²Rn decay chain) and ⁸⁵Kr. Contributions from ¹³⁶Xe double-$\beta$ decays, solar neutrino interactions, and gamma events from the materials are also taken into account. In the ROI for the WIMP search, the total ER energy spectrum is approximately flat. The distribution of ER background events in $\mathrm{cS1}$–$\mathrm{cS2}$ space is generated from the ER model fitted to ²²⁰Rn and ³⁷Ar calibration data, as discussed in Sec. II.3. The total ER rate is left unconstrained and is fitted in the final WIMP likelihood, which will be discussed in Sec. IV.

In contrast to the NR model, the uncertainty in the shape of the ER model can affect the sensitivity of the WIMP search. Ideally, the posterior distribution of all parameters should be propagated to the ER model as nuisance parameters in the WIMP search likelihood function. However, an excessive number of correlated nuisance parameters becomes computationally challenging. In XENONnT SR0, a principal component analysis (PCA) [48] is used to remove correlation among parameters from the MCMC sampler. All ER parameters shown in Tab. 2, together with $W$, $g_{1}$, and $g_{2}$ (see details in the appendix), are included in the PCA. Different from the original PCA, the variance of

$$t\equiv\sum_{i}\frac{s_{i}^{2}}{s_{i}+b_{i}}$$ (3)

along each principal component is used to quantify its importance. Here, $s_{i}$ and $b_{i}$ are the probability of a 50 GeV$/c^{2}$ WIMP and ER background in the $i^{\text{th}}$ bin of cS1–cS2 space, respectively. A larger variance of $t$ means that the uncertainty of that component can be more impactful for the WIMP sensitivity. In the end, the two most important principal components are kept as nuisance parameters. Their shape uncertainty, constrained by the calibration data, is shown in Fig. 6.

III.3 Neutron Background Model

Figure 7: Spatial (left) and cS1–cS2 (right) distribution of events tagged by the neutron veto (red markers). Black dashed and solid lines in the left plot indicate the FV and physical TPC boundary, respectively. The regions containing 68% (95%) of the total expected neutron background are indicated in dark (light) yellow. The neutron veto tagging yields four SS events (circles) and seven MS events, from which three are reconstructed with the main S2 signal inside the FV (squares) and the other four outside (diamonds). Solid and hollow markers represent the reconstructed positions of the largest and the second-largest S2 signal, respectively. The right plot only shows events inside the FV; from these, the four events within the neutron background band (large markers; one SS and three MS) were used to estimate the neutron background rate. The other three SS events are in agreement with the prediction of falsely tagged ER and surface background events (c.f. Fig. 4). For reference, the events surviving all data selections, not tagged by the neutron veto, are displayed as gray dots. Also here, the right panel displays only those data points within the FV, which correspond to the unblinded WIMP search events.

Radiogenic neutrons are mainly produced through spontaneous fission and $(\alpha,n)$ reactions in the detector materials due to intrinsic traces of radioactive impurities. The cosmogenic neutron background is subdominant [49], hence it is neglected in XENONnT SR0 analysis. The neutron yield and energy per isotope and material are calculated with the SOURCES-4A software [50], using the radioactive impurity levels of the relevant detector components obtained via the combined XENON1T and XENONnT radioassay campaigns [51, 52]. A full-chain simulation pipeline [11] is used to estimate the neutron background rate, the geometrical distribution, and the cS1–cS2 distributions for SR0.

The propagation of the neutrons in the XENONnT detector is simulated with the Geant4 toolkit [53, 54], where the recoil type and energy deposition per interaction in the target are recorded. We compute the number of photons and electrons for a given interaction via the custom-developed epix package [55], which utilizes the energy-dependent LXe response derived from the SR0 calibrations, shown in Fig. 3. This package also handles the clustering of the individual energy-depositing steps at the LXe microphysics scale before the quanta generation. The photons and electrons produced by epix are passed to the waveform simulator (WFSim) [56], which computes the S1 and S2 signals up to the waveform level, by means of a precise set of simulations- and data-driven corrections which characterize the XENONnT detector response. The event-by-event simulated waveforms share the same data structure as the science data after applying PMT and DAQ effects, which allows us to process them with the same software used for the real data (straxen [57]), as well as to apply the same data selections. The final SS neutron rate arises from weighting the rates obtained via the entire waveform simulation pipeline with the specific activities of the corresponding material and isotopic neutron yield.

When a neutron scatters in the $\sim$250 kg of LXe between the cathode and the bottom PMT array, only the scintillation light for these events can be detected. The electrons, in turn, are lost due to an electric field pointing in the opposite direction to the active volume. Neutron events that consist of scatters above and below the cathode are referred to as neutron-X events, where “X” means additional S1 contribution from charge-insensitive scatters. They are modeled as a separate background since they have a larger cS1 to cS2 ratio than normal neutron scatters.

The event parameters having discrimination power on MS interactions, such as S2 pulse shape and PMT hit patterns of S1, are matched and validated with calibrations, such that the relevant data quality selections can be applied to the simulation outputs. Notably, a validation of the multiple-to-single scatter ratio and total rate in the TPC of the AmBe calibration data was conducted, from where we obtained the systematic uncertainty associated with the accuracy of the full-chain simulations.

With this agreement between simulations and calibration data, we decided before unblinding the WIMP ROI to proceed with the sideband unblinding of the events tagged by the neutron veto. Initially, this confirmed the simulation-driven neutron background prediction. However, a mistake in the definition of the neutron veto time window was found after the unblinding of the WIMP ROI. After fixing this issue, a mismatch was found, with the neutron background rate being larger than predicted in the ROI [5]. We therefore decided to constrain the neutron background rate in a purely data-driven way based on the aforementioned sideband unblinding with the correct neutron veto tagging window. The results of the sideband unblinding are shown in Fig. 7. Based on the multiple-to-single scatter ratio of 2.2, the $53\%$ neutron veto tagging efficiency, and the three observed MS events and one SS neutron event tagged by the neutron veto in the fiducial volume, a data-driven prediction of $1.1^{+0.6}_{-0.5}$ events was derived for the SR0 exposure. We used the simulation-driven ratio between normal neutron background and neutron+X events to estimate the neutron+X event rate.

No further modification was propagated into the analysis after the data-simulations rate mismatch was identified: the $250\,\mathrm{\upmu s}$ time veto window between the TPC and the neutron veto, chosen due to the reduced background rate initially predicted, and the fiducial volume of the WIMP search remained as defined prior to unblinding. An underestimated contamination from some of the surrounding materials is considered as a possible cause of the discrepancy between the expected rate and the observed neutron background. Studies to constrain the material’s radiopurity by means of the high-energy gamma rays for specific detector regions are ongoing.

III.4 CE$\nu$NS Background Model

Neutrinos can also contribute to the NR background via coherent elastic neutrino-nucleus scattering (CE$\nu$NS). Signals from solar, atmospheric, and diffuse supernova neutrinos (DSN) will be in the WIMP search ROI. Due to the weak interaction between neutrinos and nuclei, CE$\nu$NS events are expected to be single-scatter and spatially uniform.

The recoil energy spectrum of solar CE$\nu$NS is almost identical to that of a 6 GeV$/c^{2}$ spin-independent WIMP [58], and the flux is $(5.25\pm 0.20)\times 10^{6}\,\mathrm{cm}^{-2}\mathrm{s}^{-1}$ [59, 60]. After applying all selections and their corresponding efficiencies, the expected number of events is $(0.19\pm 0.06)$ in SR0, also shown in Tab. 1. Atmospheric neutrinos and DSN mainly affect the search for heavier WIMPs. Their recoil energy spectra are taken from [49], and the SR0 expectation value is $(0.05\pm 0.02)$ events. Due to the low cross section and the similarity to WIMP interactions, CE$\nu$NS background will be the major limitation to the WIMP search sensitivity of the next generation of LXe experiments.

III.5 Accidental Coincidence Background Model

AC events consist of incorrectly paired S1 and S2 signals. These S1 and S2 signals can occur, for example, when either the S1 or the S2 signal of a physical event is not reconstructed due to detector effects, or when a single electron S2 signal is misclassified as an S1 signal. Such signals are referred to as “isolated” S1 and S2 signals. If an isolated S1 and an isolated S2 fall within the event-building time interval [8], they form an AC event.

The AC background is modeled with a data-driven approach. Isolated S1 and S2 signals as well as their surrounding S1 and S2 signals that are close in time are sampled and paired into events. S1 signals $<$150 PE are selected as isolated S1 signals, and the isolated S2 signals are taken from events whose S1 area is $<$150 PE, together with all pulses in the event window. The AC event rate is computed as

$$R_{\mathrm{AC}}=R_{\mathrm{isoS1}}\times R_{\mathrm{isoS2}}\times\Delta t,$$ (4)

where $R_{\mathrm{isoS1}}$ and $R_{\mathrm{isoS2}}$ are the isolated S1 and S2 rates, and $\Delta t$ is the event-building time interval, which is defined according to the maximum drift time of 2.2 ms in XENONnT. Occasionally during SR0, some localized high rates of S2 signals appeared in the TPC. Excluding these periods from the analysis, $R_{\mathrm{isoS1}}$ and $R_{\mathrm{isoS2}}$ remained stable throughout SR0 at an average rate of 1.5 Hz and 80 mHz, respectively. In XENONnT, the isolated S2 rate is an order of magnitude higher compared to XENON1T [61]. This can be explained by the lower electron extraction efficiency which causes an increased rate of delayed electrons. The isolated S1 and S2 signals are fed into the data processing pipeline [62, 57] to reconstruct the events. This provides the background distribution of all relevant event properties, especially the distribution in cS1 and cS2, as shown in Fig. 4 (lower right).

Validation performed on calibration data provides the rate and shape uncertainty of the AC model. For this, we compared the predicted and observed AC events in a dedicated dataset by performing both 1D and 2D GOF tests with equiprobable bins in all relevant parameters. The AC prediction is provided by the AC modeling method discussed above, and the AC dataset is selected by inverting the data selections to avoid other physical events. Because of the large statistics of ³⁷Ar calibration data, it delivers the most stringent constraint on the AC model. ³⁷Ar events with S2 areas smaller than 400 PE are selected to test the AC model. The predicted number of AC events was 731.6, while 733 events were observed in the data. The statistics in the AC model is very large so the uncertainty of the predicted AC rate is neglected in the following GOF tests. The 1D GOF tests in S1, S2, $\mathrm{Z}$, $\mathrm{R}$, and the 2D test in S1–S2 all yield p-values between 0.05 and 0.95. The result of the 2D test with a p-value of 0.61 is shown in LABEL:fig:ac_${}^{37}$Ar. The model was further validated with events removed by selections targeting the AC background inside the WIMP ROI. Similar tests were also performed on ²²⁰Rn calibration data. All these tests show the model and data to be compatible.

To suppress AC events, a selection criterion based on a gradient-boosted decision tree (GBDT) classifier utilizing the S2 pulse shape and the drift time was developed. Due to an insufficient model for the S2 pulse shape near the perpendicular wires, the GBDT classifier is only applied in the region at least 4.45 cm far from the wires (far-wire region), while in the region near the wires (near-wire region), a data-driven S2 shape selection is applied instead. The resulting AC background rate per ton-year in the near-wire region is about 5.6 times higher compared to the rate in the far-wire region. Consequently, the modeling of the TPC response is split into two parts. Due to an over-fitting issue found in the GBDT training process, we conservatively estimated the rate uncertainty to be 30% for both the near- and far-wire regions. The expectation value of AC events in the WIMP search dataset in the ROI is $(4.3\pm 0.9)$.

Figure 8: AC model validation with ³⁷Ar calibration data in S2 versus S1 area. The AC model was cross-checked against this high-statistics dataset using a 2D GOF test evaluated in an equiprobable binning scheme. The color scale shows the deviation from the predicted number of events in each bin $\mu_{\text{bin}}=14.0$ in units of sigma. The projections to S1 and S2 are shown above and to the right of the 2D plot. No significant deviations were observed in the bulk of the distribution.

III.6 Surface Background Model

Figure 9: The radial surface background model is shown for four slices in S2 area. For each slice, a Student’s generalized skew-t distribution is fitted to ²¹⁰Po $\alpha$-events. The outer radius of the physical TPC as well as of the fiducial volume are shown as a solid and a dashed black line, respectively.

In XENONnT, the TPC wall is made of PTFE, which is known to accumulate isotopes from the decay chain of atmospheric ²²²Rn, which decays down to ²¹⁰Pb, an isotope with a half-life of 22.2 years. Radioactive decays on the wall surface can result in events with reduced S2 signals due to charge losses, which gives rise to a background that can leak into the WIMP signal region. Three components contribute to the surface background in the WIMP ROI: the two $\beta$-decays of ²¹⁰Pb and ²¹⁰Bi, and the recoiling ²⁰⁶Pb following the $\alpha$-decay of ²¹⁰Po. Due to the complexity of electric field conditions near the surface and the loss of ionization electrons to the detector walls, we employ a data-driven approach to model this background component in the space of cS1–cS2–$\mathrm{R}$–$\mathrm{Z}$.

The surface background model in the space of $\mathrm{R}$–$\mathrm{Z}$ was developed using ²¹⁰Po $\alpha$-events. These events are mono-energetic (5.4 MeV) and thus easily identifiable through their characteristic S1 signal, yet they present a wide range of S2 sizes due to variable charge loss to the walls. The $\mathrm{R}$ and $\mathrm{Z}$ distributions of ²¹⁰Po $\alpha$-events were seen to match those of the lower-energy $\beta$-events. The radial profiles were modeled using a Student’s generalized skew-t distribution. The radial distribution was fitted independently for different S2 sizes as shown in Fig. 9 to account for the S2-size dependent position resolution [10]. The $\mathrm{Z}$ profile was modeled using a linear function, to account for lower rates of surface events observed at the bottom of the TPC, likely due to the increased charged insensitivity near the walls at the bottom of the detector [7].

The cS1–cS2 distribution was modeled using a 2D Gaussian adaptive kernel density estimation (aKDE), built using events reconstructed outside the detector walls. The resulting model was then validated against the events reconstructed between the walls and the fiducial volume, in order to rule out any dependence between cS1, cS2, and R. Figure 4 (lower left) shows the projection of the four-dimensional model on cS1–cS2 in the ROI. The absolute rate of events in the blinded region was inferred from the radial distributions of adjacent, non-blinded regions in cS1–cS2.

Uncertainties in the model were obtained from the $\mathrm{R}$ and $\mathrm{Z}$ fit parameter uncertainties, as those are the parameters of interest in the development of the FV. Uncertainties on the overall measured surface event rate were also propagated. Uncertainties in the (cS1, cS2) aKDE were neglected, as toy-MC tests were performed to show that they had little impact on overall expectation. For the FV, a maximum radius of $61.35\text{\,}\mathrm{cm}$ compared to the $63\text{\,}\mathrm{cm}$ used in [63] was chosen to reduce the number of surface events to $(14\pm 3)$, which improves the robustness against mismodeling.

IV Statistical Inference

In this section, we describe the statistical methods used to derive the WIMP search results, which generally follow the recommendations formulated in [46]. We first describe the likelihood function used for the search in Sec. IV.1 and its nuisance parameters in Sec. IV.2, followed by an illustration of the procedure for constructing confidence intervals and computing the discovery significance in Sec. IV.3. The GOF test to validate the best-fit models and the blinding procedure are described in Sec. IV.4. We omit some details already presented in [9], as the methods here presented build upon that previous work.

IV.1 Definition of the Likelihood Function

The fundamental ingredient for the statistical analysis of our WIMP search data is the likelihood function $\mathcal{L}(\sigma,\boldsymbol{\theta})$. It depends on the WIMP-nucleon cross section $\sigma\geq 0$, which is our parameter of interest, and a set of nuisance parameters $\boldsymbol{\theta}$. The likelihood function is defined as a product of four terms: two for the WIMP search dataset ($\mathcal{L}_{\text{near-wire}}$ and $\mathcal{L}_{\text{far-wire}}$), one for the ²²⁰Rn calibration dataset ($\mathcal{L}_{\text{cal}}$), and one for ancillary measurements of nuisance parameters ($\mathcal{L}_{\mathrm{anc}}$):

$$\mathcal{L}(\sigma,\boldsymbol{\theta})={} \mathcal{L}_{\text{near-wire}}(\sigma,\boldsymbol{\theta})$$ (5)

$$\times\mathcal{L}_{\text{far-wire}}(\sigma,\boldsymbol{\theta})$$

$$\times\mathcal{L}_{\mathrm{cal}}(\boldsymbol{\theta})$$

$$\times\mathcal{L}_{\mathrm{anc}}(\boldsymbol{\theta}).$$

In the WIMP search data, we categorize events based on their proximity to the perpendicular wires of the gate and anode electrodes, distinguishing between those that are near ($\leq$$4.45\text{\,}\mathrm{c}\mathrm{m}$, corresponding to about $17\%$ of the total FV) and far from the wires in the $\mathrm{X}$–$\mathrm{Y}$ plane. This approach allows us to account for the higher AC rate observed near the wires as discussed in Sec. III.5, without introducing a full position-dependence in the likelihood. Other backgrounds, in particular radiogenic neutron and surface background, exhibit a substantial radial dependence. Thus, the likelihood is modeled and evaluated in $\mathrm{R}$, in addition to cS1 and cS2. For the near-wire region, the radial component is not included in the modeling. Each WIMP search term is an extended unbinned likelihood function of the form

$$\mathcal{L}_{\text{region}}(\sigma,\boldsymbol{\theta})={} \mathrm{Pois}(N|\mu_{\mathrm{tot}}(\sigma,\boldsymbol{\theta}))$$ (6)

$$\times\prod_{i=1}^{N}\left[\sum_{c}\frac{\mu_{c}(\sigma,\boldsymbol{\theta})}{\mu_{\mathrm{tot}}(\sigma,\boldsymbol{\theta})}\times f_{c}(\vec{x}_{i}|\boldsymbol{\theta})\right],$$

where the index “region” runs over “near-wire” and “far-wire”. The index $i$ runs over all $N$ observed events $\vec{x}_{i}$ in the WIMP search far-wire (near-wire) dataset, where $\vec{x}_{i}$ is a vector with entries cS1, cS2, and $\mathrm{R}$ (cS1 and cS2). The PDFs $f_{c}$ of each signal and background component $c$ with expectation values $\mu_{c}$ are evaluated for each event. The total expectation value is given by $\mu_{\mathrm{tot}}(\sigma,\boldsymbol{\theta})=\sum_{c}\mu_{c}(\sigma,\boldsymbol{\theta})$.

The term $\mathcal{L}_{\mathrm{cal}}(\boldsymbol{\theta})$ is the extended unbinned likelihood function of the ²²⁰Rn calibration dataset. It depends on the two shape parameters introduced in Sec. III.2 that parameterize the range of models consistent with the ER calibration data, selected using the PCA method. By incorporating this likelihood term, we simultaneously fit the shape parameters to the ²²⁰Rn calibration dataset and the ER background model in the WIMP search dataset. Through this procedure, the constraint from the calibration data on the shape of the ER model is propagated to the final inference.

Finally, the ancillary likelihood function $\mathcal{L}_{\mathrm{anc}}(\boldsymbol{\theta})$ is a product of constraint terms for some of the nuisance parameters $\boldsymbol{\theta}$, which are detailed in the following section.

IV.2 Nuisance Parameters and Constraints

Table 1: Parameters of the XENONnT SR0 WIMP search likelihood function and their constraints. The right column shows best-fit values with an unconstrained $200\text{\,}\mathrm{GeV}\text{/}{\mathrm{\text{$c$}}}^{2}$ WIMP signal component (spin-independent coupling). For the rate parameters, all values are given in units of “events in the SR0 exposure”. Near and far-wire refers to the region in the $\mathrm{X}$–$\mathrm{Y}$ plane near and far from the perpendicular wires. In the near-wire region, which constitutes approximately $17\%$ of the total FV, the AC expectation value is comparable to the one in the far-wire region due to the higher AC rate near the wires.

Rate Parameter	Constraint	Nominal	Best Fit
ER	WIMP search data	$134.5$	$135_{-11}^{+12}$
Neutron	Ancillary measurement	$0.8\pm 0.4$	$0.8\pm 0.4$
Neutron-X	Ancillary measurement	$0.31\pm 0.16$	$0.30\pm 0.15$
CE$\nu$NS (Solar $\nu$)	Ancillary measurement	$0.19\pm 0.06$	$0.19\pm 0.06$
CE$\nu$NS (Atm + DSN)	Ancillary measurement	$0.05\pm 0.02$	$0.04\pm 0.02$
AC (near-wire)	Ancillary measurement	$2.3\pm 0.7$	$2.3\pm 0.6$
AC (far-wire)	Ancillary measurement	$2.0\pm 0.6$	$2.1\pm 0.6$
Surface	Ancillary measurement	$14\pm 3$	$12\pm 2$
220Rn calibration	220Rn dataset	$2062\pm 210$	$2052\pm 44$
Shape Parameter	Constraint	Nominal	Best Fit
ER shape parameter 1	220Rn dataset	$0.0$	$0.4\pm 0.2$
ER shape parameter 2	220Rn dataset	$0.0$	$-1.8_{-0.6}^{+0.7}$
Signal Parameter	Constraint	Nominal	Best Fit
Signal efficiency	NR model uncertainty	$1.0$	$1.0$
WIMP cross section [$10^{-47}\mathrm{cm^{-2}}$]		Parameter of interest	3.22

In addition to our parameter of interest $\sigma$, we parameterize the space of background hypotheses using nuisance parameters $\boldsymbol{\theta}$ in the likelihood function of Eq. 5. These parameters control both background rates and shape as well as the signal efficiency. In total, twelve nuisance parameters are defined, which are listed in Tab. 1. The uncertainties of the nominal expectation values correspond to the width of a Gaussian constraint term. The ancillary likelihood function $\mathcal{L}_{\mathrm{anc}}(\boldsymbol{\theta})$ is the product of all constraint terms. Most parameters are constrained via ancillary measurements, which are obtained from sidebands (e.g. data outside the ROI) or external measurements in combination with simulations. More details on the constraints on neutron, CE$\nu$NS, AC, and surface background rates were discussed in Sec. III. The ER rate is a free parameter in the fit and is fully constrained by the WIMP search data. The ER shape parameters obtained with PCA are constrained via the simultaneous fit of the ${}^{220}\mathrm{Rn}$ calibration dataset via the term $\mathcal{L}_{\mathrm{cal}}(\boldsymbol{\theta})$ in Eq. 5.

The best-fit values from the XENONnT SR0 WIMP search data [5] for a fit including an unconstrained $200\text{\,}\mathrm{GeV}\text{/}{\mathrm{\text{$c$}}}^{2}$ signal component are given in the last column of Tab. 1. The corresponding uncertainties represent the two-sided one-sigma confidence intervals derived from a profile likelihood scan in the respective parameter. More details on the construction of confidence intervals are discussed in the following.

IV.3 Confidence Intervals and Discovery Significance

The construction of confidence intervals is based on the profile likelihood test statistic

$$q(\sigma)\equiv-2\ln\frac{\mathcal{L}(\sigma,\boldsymbol{\hat{\hat{\theta}}})}{\mathcal{L}(\hat{\sigma},{\boldsymbol{\hat{\theta}}})}.$$ (7)

Quantities with a single hat denote the global maximum likelihood estimator of a parameter, while $\boldsymbol{\hat{\hat{\theta}}}$ denotes the set of nuisance parameters that maximize the likelihood if the cross section is fixed to a given $\sigma$. The likelihood is defined for a specific signal model, for example a $\mathrm{WIMP}$ with a certain mass and interaction type (e.g. spin-independent or spin-dependent coupling). In our statistical inference, we consider signal models across a range of WIMP masses, from $6\text{\,}\mathrm{GeV}\text{/}{\mathrm{\text{$c$}}}^{2}$ to $500\text{\,}\mathrm{GeV}\text{/}{\mathrm{\text{$c$}}}^{2}$. For higher WIMP masses, the PDF remains constant and the flux for a given cross section decreases approximately linearly with mass, as discussed in Sec. III.1. According to this, limits can be extrapolated even beyond $500\text{\,}\mathrm{GeV}\text{/}{\mathrm{\text{$c$}}}^{2}$.

Knowing the distribution of $q(\sigma)$ under different hypotheses is essential for calculating discovery significances and confidence intervals. Due to the low rate of background resembling the signal, asymptotic formulae as discussed in [64] are not applicable. Instead, toy-MCs are used to estimate the distribution $g(q(\sigma)|\sigma)$ of the test statistic given that $\sigma$ is the true cross section. We use the custom-developed blueice [65] framework to define the likelihood function. This Python package provides an efficient interpolation (“template morphing”) between PDFs evaluated for different discrete values of shape nuisance parameters. This allows shape parameters to be considered for which we have no analytical description.

For each signal model considered, we compute the discovery significance and the confidence interval by comparing the measured test statistic $q(\sigma)$ with the distribution under each hypothesis $g(q(\sigma)|\sigma)$ [66, 67]. Testing $\sigma=0$ yields the (local) discovery significance, while confidence intervals are computed by finding the region where $\sigma$ is rejected at a $90\%$ confidence level (CL) given the data, which is illustrated for three background-only toy-MCs (three black parabola-like lines) in Fig. 10 (top). The minimum of each curve corresponds to the respective best-fit value $\hat{\sigma}$. Upper limits (ULs) and lower limits are obtained by finding the cross sections at which $q(\sigma)$ reaches the critical region. The threshold of the 90 $\%$ CL critical region is precomputed as the 90th percentile of the test statistic distribution from toy-MCs with injected signal corresponding to the respective cross section. For low cross sections, it deviates from the asymptotic two-sided threshold indicated as a horizontal gray dotted line.

Repeating this procedure $\mathcal{O}(1000)$ times yields the expected distribution of ULs that can be obtained in the absence of any signal, shown in Fig. 10 (bottom). The experiment’s sensitivity is given by the median of these ULs and the sensitivity band (“Brazil band”) marked by the $\pm$2-sigma and $\pm$1-sigma quantiles indicates the spread. The distribution of lower limits for background-only toy-MCs is also shown in Fig. 10 (bottom). As recommended in [46], we decide on a discovery significance threshold for reporting two-sided confidence intervals of $3$ sigma, following [9].

Figure 10: Illustration of the confidence interval construction and distribution of limits on the WIMP-nucleon cross section $\sigma$ of a $200\text{\,}\mathrm{GeV}\text{/}{\mathrm{\text{$c$}}}^{2}$ WIMP for background-only toy-MCs. Top: The test statistic $q(\sigma)$ as a function of the cross section $\sigma$ is shown for three toy-MCs. The intersections with the threshold of the critical region (gray line) yield the 90 $\%$ CL upper (blue diamonds) and lower limits (purple diamonds). Bottom: (Inverse) cumulative distribution function (CDF) of upper (lower) limits for background-only toy-MCs. The bands containing 68 $\%$ (green) and 95 $\%$ (yellow) of ULs, as well as the median UL (dashed line), are indicated. Note the shared abscissa with the cross section (bottom) and the corresponding expected number of signal events (top).

Statistical fluctuations as well as mismodeling, such as overestimated background rates, can yield arbitrarily low ULs, which may result in the spurious exclusion of models beyond the experiment’s sensitivity. To mitigate this issue, various methods have been proposed [68, 69, 70]. The XENONnT SR0 WIMP search follows the recommendation of [46] to use power-constrained limits (PCL) [70]. In this method, a signal size threshold is selected, at which the experiment has a “rejection power” $M_{\mathrm{min}}$ – the probability of excluding a given signal under the background-only hypothesis. If an UL falls below this threshold, the threshold value is reported instead. These thresholds correspond to the quantiles of the UL distribution used to compute the sensitivity band illustrated in Fig. 10. For instance, choosing $M_{\mathrm{min}}=0.16$ sets the threshold to the -1 sigma quantile of the band, while $M_{\mathrm{min}}=0.5$ truncates ULs at the band’s median. In [46], the fiducial choice of $M_{\mathrm{min}}=0.16$ was suggested. However, the choice was erroneously based on the discovery power, which corresponds to the probability of rejecting the background-only hypothesis given an alternative hypothesis with a specific signal size. In the absence of an updated recommendation, we chose the very conservative choice of $M_{\mathrm{min}}=0.5$ for the first WIMP search of XENONnT.

IV.4 Goodness of Fit and Blinding Procedure

To verify that the best-fit models describe our WIMP search data well, we defined a GOF test before unblinding, and after studies to optimize the power to reject mismodeling. The test uses a binned Poisson likelihood $\chi^{2}$ with 15 equiprobable bins in the cS1–cS2 space, defined from the model being tested. To compute a p-value, the distribution of the $\chi^{2}$ test statistic under the best-fit hypothesis is derived through toy-MCs. Specifying the test and acceptance threshold (90% CL) before unblinding ensured that the results were not influenced by statistical fluctuations expected from the low-statistic WIMP search dataset. For the background-only fit we found a p-value of $0.67$ and for a fit with an additional $200\text{\,}\mathrm{GeV}\text{/}{\mathrm{\text{$c$}}}^{2}$ WIMP signal the p-value is $0.63$. Both indicate no strong mismodeling in the predefined parameter space. Using GOF tests with well-studied power to discover signal-like mismodeling was also designed to replace the “safeguard” ER shape parameter defined in [71], due to its computational cost and susceptibility to some kinds of mismodeling.

When searching for new signals in data, it is crucial to avoid the experimenter’s bias on the result [72]. In our WIMP search, we adopt a common strategy for bias mitigation: blinding the signal region until all steps of the data analysis are finalized. Initially, events in the ER and NR bands with a reconstructed ER energy below $20\,\mathrm{keV_{ER}}$ were blinded. This involved defining the ER and NR bands in cS1 and cS2 based on quantiles in cS2 for slices in cS1 of calibration data. In the first step, all events above $10\,\mathrm{keV_{ER}}$ and those above the $-2$ sigma quantile of the ER band were unblinded for the analysis presented in [63]. The signal region for the WIMP search (indicated in Fig. 5) remained blinded until the analysis procedure was finalized. The unblinded SR0 WIMP search data showed no significant excess for any of the tested WIMP masses with local discovery p-values $\geq 0.2$. For this reason, we reported new ULs on the spin-independent WIMP-nucleon cross section across WIMP masses ranging from $6\text{\,}\mathrm{GeV}\text{/}{\mathrm{\text{$c$}}}^{2}$ to $500\text{\,}\mathrm{GeV}\text{/}{\mathrm{\text{$c$}}}^{2}$, with a minimum of $2.58\times 10^{-47}\;\mathrm{cm}^{2}$ at $28\text{\,}\mathrm{GeV}\text{/}{\mathrm{\text{$c$}}}^{2}$.

V Summary

The XENONnT WIMP dark matter search relies on a detector response model as well as simulation- and data-driven background models. These were combined to construct a statistical model in cS1, cS2, and $\mathrm{R}$, which was used to infer limits on WIMP-nucleus scattering cross sections.

The full detector response model for electronic recoil (ER) and nuclear recoil (NR) interactions, including both the xenon emission and detector reconstruction model, was successfully fitted to calibration data. Accurate simulations of particle interactions up to the data acquisition waveform level made this possible, in particular, to correctly model the S2 multiplicity of events with several, potentially unresolved energy deposits. Using these models, we derived the distributions for ER and NR backgrounds, as well as the signals, in our analysis space. Except for ER, the background rates were constrained with ancillary measurements. The radiogenic neutron background rate was constrained by first matching the simulated ratio of multiple-to-single scatter interactions and the neutron veto tagging efficiency with NR calibration data. After unblinding the neutron veto-tagged events, these three inputs were combined to derive a prediction for the remaining non-vetoed single-scatter neutron background. Two shape parameters were propagated to the final inference to account for the uncertainty of the ER model in cS1–cS2 space. Accidental coincidence and surface background were modeled with data-driven approaches. The validity of the models was confirmed in calibration data as well as science data outside the region of interest of the WIMP search.

A blind analysis was performed for the first science data from XENONnT. The statistical methods largely follow the previous XENON1T approach and community recommendations, by using a toy-MC-calibrated profile log-likelihood ratio test statistic. One departure from these recommendations was raising the minimal required power of the power-constrained limit (PCL) threshold from 0.15 to 0.5, corresponding to placing upper limits only at or above the median unconstrained upper limit. Both the calibration fits and the final best-fit model were assessed and found acceptable with goodness-of-fit tests that were chosen based on their mismodeling rejection power, and defined prior to unblinding the data. Analysis of the data with an exposure of 1.1 tonne-years revealed no signal excess over backgrounds. Therefore, new upper limits on the spin-independent WIMP-nucleon cross section were derived.

Appendix A Details on the ER Emission Model

Let $E$ be the recoil energy of an ER event. The total number of detectable quanta $N_{\mathrm{q}}$ is sampled from a normal distribution

$$N_{\mathrm{q}}\sim\mathrm{Norm}(E/W,\sqrt{fE/W}).$$ (8)

$W$ is the mean energy to generate a quantum and $f$ is the Fano factor $\leavevmode\nobreak\ 0.059$ [73, 74]. The ER energy deposit will produce both excitons and electron-ion pairs. Using the mean ratio between number of excitons and ions $\left\langle N_{\mathrm{ex}}/N_{\mathrm{i}}\right\rangle$, we simulate their numbers from a binomial distribution:

$$N_{\mathrm{i}}\sim\mathrm{Binom}\left(N_{\mathrm{q}},\frac{1}{1+\langle N_{\mathrm{ex}}/N_{\mathrm{i}}\rangle}\right),$$ (9)

$$N_{\mathrm{ex}}=N_{\mathrm{q}}-N_{\mathrm{i}}.$$ (10)

A fraction $r$ of ions recombine with electrons, depending on the electric drift field $F$. We parametrize the mean value of $r$ in the same way as in XENON1T [9]:

$$\langle r\rangle=\frac{1}{\mathrm{e}^{-(E-q_{0})/q_{1}}+1}\left(1-\frac{\log(1+\langle N_{\mathrm{i}}\rangle\varsigma)}{\langle N_{\mathrm{i}}\rangle\varsigma}\right),$$ (11)

where

$$\varsigma=\frac{1}{4}\gamma\mathrm{e}^{-E/\omega}F^{-\delta}.$$ (12)

The fluctuation of $r$ is parametrized via

$$\Delta r=q_{2}(1-\mathrm{e}^{-E/q_{3}}),$$ (13)

and the true recombination fraction $r$ is then sampled from

$$r\sim\mathrm{Norm}(\langle r\rangle,\Delta r).$$ (14)

All fitted ER parameters are listed in Tab. 2. Finally, the numbers of produced electrons and photons are given by

$$N_{\mathrm{e}}\sim\mathrm{Binom}(N_{\mathrm{i}},1-r),$$ (15)

$$N_{\gamma}=N_{\mathrm{q}}-N_{\mathrm{e}}.$$ (16)

Appendix B Details on the NR Emission Model

Let $E$ be the recoil energy of an NR event. The total number of produced quanta $N$ is

$$N\sim\mathrm{Norm}(E/W,\sqrt{fE/W}).$$ (17)

In contrast to electron recoils, recoiling xenon nuclei lose a significant amount of their recoil energy via atomic motion and collisions with other xenon atoms, which are processes that are undetectable in a LXe TPC. This quanta loss is modeled following the Lindhard theory [75], with the so-called Lindhard quenching factor $L$, such that the number of detectable quanta becomes

$$N_{\mathrm{q}}\sim\mathrm{Binom}(N,L).$$ (18)

Following the parametrization in [14], the Lindhard factor is given by

$$L=\frac{\kappa\,g(\epsilon)}{1+\kappa\,g(\epsilon)},$$ (19)

where $g(\epsilon)$ is a function of deposited energy via

$$g(\epsilon)=3\epsilon^{0.15}+0.7\epsilon^{0.6}+\epsilon,$$ (20)

$$\epsilon=11.5\,Z^{-7/3}(E/\mathrm{keV}),$$ (21)

with the nuclear charge number $Z=54$ for xenon. The numbers of produced ions and excitons are then simulated by

$$N_{\mathrm{i}}\sim\mathrm{Binom}\left(N_{\mathrm{q}},\frac{1}{1+\langle N_{\mathrm{ex}}/N_{\mathrm{i}}\rangle}\right),$$ (22)

$$N_{\mathrm{ex}}=N_{\mathrm{q}}-N_{\mathrm{i}},$$ (23)

with the exciton-to-ion ratio parametrized as

$$\langle N_{\mathrm{ex}}/N_{\mathrm{i}}\rangle=\alpha F^{-\zeta}(1-\mathrm{e}^{-\beta\epsilon}).$$ (24)

Unlike ERs, the recombination fluctuation in NRs is usually small, thus the number of photons produced from the recombination of electron-ion pairs is

$$N_{\gamma}^{\mathrm{re}}\sim\mathrm{Binom}(N_{\mathrm{i}},r),$$ (25)

where $r$ follows the Thomas-Imel box model [76]

$$r=1-\frac{\log(1+N_{i}\varsigma)}{N_{i}\varsigma},$$ (26)

$$\varsigma=\gamma F^{-\delta}.$$ (27)

The number of excitons is further reduced by bi-excitonic and Penning quenching effects [77], such that the number of photons produced from de-excitation becomes

$$N_{\gamma}^{\mathrm{de}}\sim\mathrm{Binom}(N_{\mathrm{ex}},f_{l}),$$ (28)

with the scintillation quenching factor

$$f_{l}=\frac{1}{1+\eta\epsilon^{\lambda}}.$$ (29)

Then the final numbers of photons and electrons are

$$N_{\gamma}=N_{\gamma}^{\mathrm{de}}+N_{\gamma}^{\mathrm{re}},$$ (30)

$$N_{\mathrm{e}}=N_{\mathrm{i}}-N_{\gamma}^{\mathrm{re}}.$$ (31)

Appendix C Details on the Detector Response Model

The simulations of $\mathrm{S1}$ and $\mathrm{S2}$ signals from $N_{\gamma}$ and $N_{\mathrm{e}}$ are almost independent of each other.

To simulate $\mathrm{S1}$ signals, we first introduce the spatially dependent scintillation gain $\tilde{g_{1}}$ with DPE excluded:

$$\tilde{g_{1}}(\mathrm{X},\mathrm{Y},\mathrm{Z})=\frac{g_{1}}{1+p_{\mathrm{DPE}}}\cdot\mathrm{S1corr}^{-1}(\mathrm{X},\mathrm{Y},\mathrm{Z}),$$ (32)

where $g_{1}$ is the averaged scintillation gain, $\mathrm{S1corr}$ is the relative spatial correction factor, and $p_{\mathrm{DPE}}$ is the probability of the double photoelectron emission effect [17]. Then, for a given number of photons $N_{\gamma}$ and the position of the recoil $\mathrm{X}$, $\mathrm{Y}$, $\mathrm{Z}$, the number of photons detected by the PMTs is

$$N_{\mathrm{PhD}}\sim\mathrm{Binom}(N_{\gamma},\tilde{g_{1}}(\mathrm{X},\mathrm{Y},\mathrm{Z})).$$ (33)

Now accounting for the DPE effect, the expected number of photoelectrons (PE) in the S1 signal is simulated via

$$N_{\mathrm{DPE}}\sim\mathrm{Binom}(N_{\mathrm{PhD}},p_{\mathrm{DPE}}),$$ (34)

$$N_{\mathrm{S1,PE}}\sim N_{\mathrm{PhD}}+N_{\mathrm{DPE}}.$$ (35)

However, the $\mathrm{S1}$ is not always equal to the number of PEs due to reconstruction bias. This is modeled as

$$\mathrm{S1}/N_{\mathrm{S1,PE}}-1\sim\mathrm{Norm}(\delta_{\mathrm{S1}},\Delta_{\mathrm{S1}}).$$ (36)

Here, both $\delta_{\mathrm{S1}}$, $\Delta_{\mathrm{S1}}$ are functions of $N_{\mathrm{PhD}}$ obtained from the XENONnT MC. In the end, the relative spatial correction is applied back to the $\mathrm{S1}$ signal, yielding

$$\mathrm{cS1}=\mathrm{S1}\cdot\mathrm{S1corr}(\mathrm{X}^{\prime},\mathrm{Y}^{\prime},\mathrm{Z}),$$ (37)

where $\mathrm{X}^{\prime},\mathrm{Y}^{\prime}$ are the event positions as reconstructed in data, which are smeared by the S2-size-dependent resolution of the position reconstruction.

To get the $\mathrm{S2}$ signal from $N_{\mathrm{e}}$, the first step is to simulate the electron loss during the drift. The survival probability due to attachment to electronegative impurities in LXe is given by

$$p_{\mathrm{loss}}(\mathrm{Z})=\mathrm{e}^{-\mathrm{Z}/(\tau v)},$$ (38)

where $\tau$ is the electron drift survival time (“electron lifetime”) and $v$ is the drift velocity of electrons in LXe. Inside the FV the electric field is uniform enough to approximate $v$ as a constant. Electron losses due to drift field effects close to the wall are accounted for via a spatially dependent charge-insensitive-volume (CIV) probability function $p_{\mathrm{CIV}}(\mathrm{R},\mathrm{Z})$ from field simulations of XENONnT [7]. The number of surviving electrons is then given by

$$N_{\mathrm{surv}}\sim\mathrm{Binom}\left(N_{\mathrm{e}},\>p_{\mathrm{loss}}(\mathrm{Z})\cdot p_{\mathrm{CIV}}(\mathrm{R},\mathrm{Z})\right).$$ (39)

At the liquid-gas interface, a fraction of electrons are extracted,

$$N_{\mathrm{extr}}\sim\mathrm{Binom}(N_{\mathrm{surv}},\epsilon(\mathrm{X},\mathrm{Y})).$$ (40)

The extraction efficiency is $\mathrm{X}$-$\mathrm{Y}$ dependent and can be calculated by

$$\epsilon(\mathrm{X},\mathrm{Y})=g_{2}\cdot\mathrm{S2corr}^{-1}(\mathrm{X},\mathrm{Y})/G,$$ (41)

where $g_{2}$ is the averaged ionization gain, $\mathrm{S2corr}$ is the relative spatial correction factor, and $G\sim 31\,\mathrm{PE/e^{-}}$ is the single electron gain. Assuming the fluctuation of the secondary scintillation process is Poisson-like, the expected number of $\mathrm{S2}$ PEs is

$$N_{\mathrm{S2,PE}}\sim\mathrm{Norm}(N_{\mathrm{extr}}G,\sqrt{N_{\mathrm{extr}}G}).$$ (42)

Similar to $\mathrm{S1}$, the $\mathrm{S2}$ reconstruction bias is modeled by

$$\mathrm{S2}/N_{\mathrm{S2,PE}}-1\sim\mathrm{Norm}(\delta_{\mathrm{S2}},\Delta_{\mathrm{S2}}).$$ (43)

Finally, the correction is applied to the $\mathrm{S2}$,

$$\mathrm{cS2}=\mathrm{S2}\cdot\mathrm{S2corr}(\mathrm{X}^{\prime},\mathrm{Y}^{\prime})\cdot\mathrm{e}^{\mathrm{Z}/\tau^{\prime}v}.$$ (44)

Here, $\tau^{\prime}$ is the electron lifetime used in the correction, which could be slightly off from the true value. The difference is one of the parameters to fit in the band fitting.

For both S1 and S2 signals, data quality selection efficiencies are applied which depend on the respective signal sizes. For S1 signals, we additionally apply the S1 reconstruction efficiency which is a function of the number of detected photons $N_{\mathrm{PhD}}$.

For NR multi-scatter events, the simulation is slightly altered. Here the total number of detected photons is defined as the sum of the detected photons for each single energy deposition $i$,

$$N_{\mathrm{PhD}}=\sum_{i}N_{\mathrm{PhD}}^{(i)}.$$ (45)

Similarly, we sum over the number of S2 PEs of only those scatters that contribute to the S2,

$$N_{\mathrm{S2,PE}}=\sum_{i_{\mathrm{S2}}}N_{\mathrm{S2,PE}}^{(i_{\mathrm{S2}})}.$$ (46)

Reconstruction biases, signal corrections, and selection efficiencies are correspondingly applied to the merged signals.

Appendix D Details on the Calibration Fit Results

The parameters in the ER and NR band fits are summarized in Tab. 2, showing both the applied prior and the marginal posterior values. Parameters where prior values are given without uncertainties are fixed in the fit. The prior for $\langle N_{\mathrm{ex}}/N_{\mathrm{i}}\rangle$ is a uniform distribution in the stated range. For the posteriors, the stated values and uncertainties correspond to the median and the central 68% of the marginal posterior distributions.

Table 2: Prior and marginal posterior distributions of each parameter in the ER and NR emission models.

Parameter Prior Marginal posterior Unit $W$ $13.7\pm 0.2$ $13.7^{+0.2}_{-0.2}$ eV $f$ 0.059 0.059 - ER parameters $\langle N_{\mathrm{ex}}/N_{\mathrm{i}}\rangle$ $0.06$ - $0.20$ $0.13^{+0.04}_{-0.04}$ - $\gamma$ free $0.13^{+0.03}_{-0.02}$ - $\delta$ free $0.34^{+0.07}_{-0.07}$ - $\omega$ free $57^{+15}_{-12}$ keV $q_{0}$ free $1.32^{+0.17}_{-0.20}$ keV $q_{1}$ free $0.47^{+0.07}_{-0.05}$ keV $q_{2}$ free $0.030^{+0.002}_{-0.002}$ - $q_{3}$ free $0.47^{+0.40}_{-0.31}$ keV NR parameters $\zeta$ 0.047 0.047 - $\delta$ 0.062 0.062 - $\alpha$ free $0.92^{+0.07}_{-0.06}$ - $\gamma$ free $0.016^{+0.001}_{-0.001}$ - $\beta$ $239^{+56}_{-18}$ $334^{+40}_{-43}$ - $\kappa$ $0.139^{+0.006}_{-0.005}$ $0.138^{+0.005}_{-0.006}$ - $\eta$ $3.3^{+10.6}_{-1.4}$ $10.0^{+6.8}_{-5.9}$ - $\lambda$ $1.14^{+0.90}_{-0.18}$ $1.40^{+0.61}_{-0.38}$ -