First Constraints on WIMP-Nucleon Effective Field Theory Couplings in an Extended Energy Region From LUX-ZEPLIN

Following the first science results of the LUX-ZEPLIN (LZ) experiment, a dual-phase xenon time projection chamber operating from the Sanford Underground Research Facility in Lead, South Dakota, USA, we report the initial limits on a model-independent non-relativistic effective field theory describing the complete set of possible interactions of a weakly interacting massive particle (WIMP) with a nucleon. These results utilize the same 5.5 t fiducial mass and 60 live days of exposure collected for the LZ spin-independent and spin-dependent analyses while extending the upper limit of the energy region of interest by a factor of 7.5 to 270 keV nr . No significant excess in this high energy region is observed. Using a profile-likelihood ratio analysis, we report 90% confidence level exclusion limits on the coupling of each individual non-relativistic WIMP-nucleon operator for both elastic and inelastic interactions in the isoscalar and isovector bases.

Following the first science results of the LUX-ZEPLIN (LZ) experiment, a dual-phase xenon time projection chamber operating from the Sanford Underground Research Facility in Lead, South Dakota, USA, we report the initial limits on a model-independent non-relativistic effective field theory describing the complete set of possible interactions of a weakly interacting massive particle (WIMP) with a nucleon.These results utilize the same 5.5 t fiducial mass and 60 live days of exposure collected for the LZ spin-independent and spin-dependent analyses while extending the upper limit of the energy region of interest by a factor of 7.5 to 270 keVnr.No significant excess in this high energy region is observed.Using a profile-likelihood ratio analysis, we report 90% confidence level exclusion limits on the coupling of each individual non-relativistic WIMP-nucleon operator for both elastic and inelastic interactions in the isoscalar and isovector bases.

I. INTRODUCTION
Current and next-generation dark matter (DM) direct detection experiments searching for weakly interacting massive particles (WIMPs) will reach unprecedented levels of exposure during their operational lifetimes, allowing them to probe new parameter space for WIMPs, as demonstrated by the first LZ [1] and XENONnT [2] spin-independent (SI) and spin-dependent (SD) WIMPnucleon interaction limits.Both the SI and SD interaction models assume that the momentum dependence of the interaction is suppressed due to the non-relativistic velocity of the WIMP in the nucleus rest frame, hence allowing for the use of a zero momentum cross-section [3].However, this may be an oversimplification of the dynamics of the interaction as, despite being small, the velocity is still non-zero.By considering interactions with momentum dependence, it is possible to increase the potential for the discovery of WIMP-nucleon interactions by direct detection experiments.
The rate of a more generalized interaction may depend on the momentum transfer; if the momentumindependent components are suppressed, as in the interaction of a nucleon and DM with a magnetic dipole moment, the dominant component will be momentum-dependent.When considering all possible interactions in the nucleon frame, the momentum may be non-negligible and is certainly non-negligible for WIMP-parton interactions [4].Given the variety of possible WIMPnucleon interactions, the most comprehensive approach is to consider the WIMP-nucleon interaction in a modelindependent way.Since the non-gravitational interactions of DM are unknown, it is compelling to describe the interaction through an effective field theory (EFT) that captures a significant amount of the possible physics [5].Constraints on a comprehensive set of WIMP-nucleon interactions generated using an EFT [6,7], have previously been established by PandaX-II [8], LUX [9], XENON100 [10], DEAP-3600 [11], Darkside-50 [12] and PICO-60 [13].This report further constrains the coefficients within this EFT by applying a profile likelihood ratio analysis to the first science data of the LZ experiment.An extended energy window is used to increase the coverage of the signal model profiles, which may produce significant interaction rates at energies beyond that of a typical SI or SD WIMP search.As such, the calibrations, efficiencies, and relevant backgrounds are reassessed for this wider energy range.

A. Non-Relativistic Effective Field Theory
A full derivation of the WIMP-nucleon non-relativistic effective field theory (NREFT), which describes all possible interactions between WIMPs and a target nucleus, has been developed by Fan et al. [6] and Fitzpatrick et al. [7].In an NREFT, one expands the effective Lagrangian in powers of the momentum transferred in the WIMP-nucleus interaction.The cutoff scale for this theory is dictated by the maximum momentum transfer, around 200 MeV, set by the reduced mass of the WIMPnucleus system and escape velocity of the WIMP [14].The pion degrees of freedom are integrated out of this EFT, restricting its validity to energies below the mass of the pion [7].While this NREFT reliably describes all possible interactions at these energies, in order to match it to an ultraviolet theory of dark matter, the QCD dynamics that were ignored during the creation of the NREFT have to be brought back.Examples of doing this are described in Refs.[15][16][17] and experimental limits on WIMP-pion interactions are described in Ref. [18].
In the Fitzpatrick et al.NREFT [7], the corresponding complete interaction Lagrangian can be written in terms of linear combinations of effective operators, each multiplied by a coefficient.These operators can be expressed in terms of WIMP-nucleon couplings, which are proportional to low-energy constants (parameters such as the WIMP axial-vector or scalar coupling constant) in both the effective Lagrangian and the nuclear form factor.The nuclear form factor describes the spatial distribution of the nuclear density and is calculated using models that consider the structure and dynamics of the target nucleus.The NREFT is thus a powerful tool for interpreting experimental data in a model-independent way.Furthermore, the power counting [19] scheme of NREFT allows for a systematic and controlled expansion of the scattering cross-section in terms of the momentum transfer, which is essential for calculating the sensitivity of direct detection experiments to WIMP dark matter.The derivation of the NREFT introduces both momentumdependent and velocity-dependent operators to describe all possible mechanisms for the WIMP-nucleon interaction.By enforcing momentum conservation and Galilean invariance [7], the operators can be reduced to a basis of four Hermitian quantities: where ⃗ q is the momentum transferred from the WIMP to the nucleon, m N is the nucleon mass, ⃗ v ⊥ is the component of the relative velocity between the WIMP and the nucleon that is perpendicular to the momentum transfer, ⃗ S χ is the spin of the WIMP, and ⃗ S N is the spin of the relevant nucleon.By forming linear combinations of the Hermitian quantities, up to second-order in ⃗ q, a total of fifteen independent and dimensionless NREFT operators can be constructed [14]: (2) As only non-relativistic interactions are considered in this analysis, any operator with a quadratic or greater dependency on v ⊥ is not included (namely O 2 ).This analysis is conducted in the isoscalar (s) and isovector (v) bases as opposed to the proton and neutron bases used in some previous analyses.It can be argued that proton and neutron are much more natural bases to produce WIMP-nucleon scattering constraints, as this is what the WIMP interacts with.However, by exploiting the fact that isospin is a good quantum number in nuclear systems, it is possible to decompose the fundamental constraints on the operators when comparing experiments with varying target nuclei.An additional benefit to using the isoscalar and isovector bases is that it is possible to compare to the limits set on the elastic and inelastic WIMP-nucleon NREFT operators by the XENON collaboration [10], as well as the inelastic limits set by LUX [9].
This report also considers inelastic scatters where the masses of the incoming and outgoing DM particles can differ, allowing for the WIMP to transition into a more massive state during the scattering interaction [20].This mass difference, denoted by δ m ≡ m χ,out − m χ,in , results in a minimum required recoil energy to produce an inelastic interaction.This in turn leads to a suppressed recoil rate at lower energies and an observed signal that is concentrated at higher energies.In some cases, where elastic scattering is suppressed, inelastic transitions between WIMP states become the primary method of interaction [21,22].A slight modification of the Hermitian basis vectors is required to study the case of inelastic WIMP-nucleon interactions.In elastic interactions, energy conservation results in ⃗ v ⊥ • ⃗ q = 0. To preserve energy conservation in inelastic interactions where there is a nonzero mass splitting, δ m , the following condition must be satisfied [23]: where µ N is the reduced mass of the WIMP-nucleon system [23].This requirement is incorporated into the Hermitian basis vectors by replacing the perpendicular ve-locity in Eq. 1 with where v ⊥ is the perpendicular velocity for elastic scattering.A similar replacement is made for each operator, O i .This report considers mass splitting values in the 0-250 keV range, which are well-motivated in many WIMP models [20,23].
Before experimental limits can be translated into bounds on the operator coefficients of Eq. 2, it is necessary to calculate nuclear response functions for DM elastic scattering.It is common practice to use the fullbasis shell-model calculations of each contributing Xe isotope using the GCN5082 interaction when determining these responses [24].Recently, these calculations have been revised, leading to a much-improved determination of the one-body nuclear density matrices [25].In this analysis, the updated density matrices were incorporated into DMFormfactor-v6 [14], a software package for calculating the NREFT dark matter-nucleus scattering.These response functions were subsequently integrated into WimPyDD, a framework for modeling the interaction of dark matter with atomic nuclei in direct detection experiments [26].This modified version of WimPyDD was validated against both DMFormfactor-v6 and the published interaction Lagrangian recoil spectra by Pan-daX and Haxton [8].Fig. 1 shows the differential rate spectra, generated using the modified WimPyDD, for each non-relativistic operator considered in this analysis, assuming a WIMP mass of 1000 GeV/c 2 .

II. DETECTOR, DATA AND SELECTION
A. The LZ Detector The LZ experiment is housed in the Davis Campus of the Sanford Underground Research Facility, 4850 (4300 m.w.e) below Lead, South Dakota, USA, which provides an attenuation of the cosmic muon flux by a factor of 10 6 [27,28].LZ utilizes a low-background, dual-phase time projection chamber (TPC) which contains 7 tonnes of liquid xenon (LXe) in the active volume [27,28].The cylindrical TPC, with both a height and diameter of approximately 1.5 m, is equipped with two arrays consisting of 251 and 243 R11410-22 photomultiplier tubes (PMTs) at the top and bottom of the detector, respectively, to detect scintillation light.
Energy depositions in the LXe volume are detected by collection of the scintillation produced from the initial interaction (S1), as well as from the electroluminescent light from extracted electrons (S2).Electrons are collected due to the presence of electric fields applied across the detector volume: a 192 V/cm drift field and a 7.3 kV/cm gas electroluminescence field.The S1 and S2 signals are reported in units of photons detected (phd), accounting for the double photon emission effect, which can occur when VUV light is incident on a PMT.A requirement that scintillation light is detected in at least three PMTs (three-fold coincidence) is used when defining an S1 in the data.
Additionally, the S2:S1 ratio differs for electronic recoils (ERs) and nuclear recoils (NRs), providing discrimination power between different interaction types.This is illustrated in Fig. 2 using LZ calibration datasets.The position-corrected S1 and S2 observables are typically used (S1c and S2c, respectively), to improve the energy resolution and discrimination between ERs and NRs.Position-corrected observables also allow for the use of linear scaling factors, g 1 and g 2 , to correlate the S1c and S2c signals, respectively, to the original number of photons and electrons produced (n ph and n e ): To reject internal and external backgrounds, the LZ experiment includes two veto detectors: a xenon "Skin" veto surrounding the active mass, which is outfitted with 93 1-inch and 38 2-inch PMTs; and a near-hermetic "outer detector" (OD) consisting of acrylic tanks containing 17 tonnes of gadolinium-loaded liquid scintillator (0.1% by mass).The skin is designed to identify multiple scattering interactions entering or exiting the TPC, while the outer detector is designed to capture and identify neutrons that may scatter in the TPC.The entire system is housed in a tank that is filled with 238 tonnes of ultra-pure water, which provides shielding from ambient radioactive backgrounds emitted by the cavern rock [27,28].120 8-inch PMTs are situated around the walls of this tank to observe any light produced in the OD.
A 60 live-day exposure using a 5.5 ± 0.2 tonne fiducial volume (FV) was collected between December 2021 and May 2022.Using events from a dataset with S1c between 3 and 80 phd, uncorrected S2 > 600 phd and log 10 (S2c) < 5, LZ set world-leading constraints on SI WIMP-nucleon interactions for masses greater than 9 GeV/c 2 , with the most stringent limit excluding crosssections greater than 9.2 × 10 −48 cm 2 for 36 GeV/c 2 WIMPs [1].This report extends the region of interest (ROI) of the same dataset to include S1c signals between 3 and 600 phd to give sensitivity to NREFT interactions with their most significant rate contribution outside the typical search window used when considering SI and SD interactions.The maximum S2c considered in this analysis is set at log 10 (S2c) = 4.5 to remove ER backgrounds far from the NR signal region, given that leakage from ERs becomes less significant at higher energies (see Fig. 2).The lower bound on uncorrected S2 is maintained at 600 phd.

′′
(bottom-center) and all other responses (bottom-right).The spectra were generated with a coupling strength of unity, excluding the possibility of interference terms.The shaded gray regions indicate the energies at which the detection efficiency is below 50% after all data analysis cuts have been applied (see Fig. 4).

B. Calibrations and Data Selections
As described in Ref. [1], the ER and NR responses were measured using dedicated in-situ calibrations with tritiated methane (0-18.6 keV ERs) and D-D fusion neutrons (0-74 keV NRs).The detector and model response parameters from NEST 2.3.7 (Noble Element Simulation Technique) [29,30] were tuned to the ER and NR calibration data in order to reproduce the observed data.
This tuning was used along with constraints from the energy reconstruction of several monoenergetic peaks from background and calibration sources to find the S1c and S2c scaling factors defined in Eq. 5: g 1 = 0.114 ± 0.002 phd/photon and g 2 = 47.1 ± 1.1 phd/electron.
Extending the ER and NR response models through the extended energy ROI is one of the key challenges in performing an NREFT analysis.For this, β emissions from 212 Pb following an injection of 220 Rn provide a calibration of the ER response through the 3 -600 phd NREFT search window.Reproducing the 212 Pb response does not require altering g 1 or g 2 from Ref [1]; however, NEST underestimates the observed ER band widths above 100 phd, worsening at higher energies (O(10%) disagreement at 600 phd).To account for this, the functionality in NEST for energy-dependent smearing of pulse areas is utilized, similar to the methods reported in Ref. [31].This results in proper reproduction of the ER response in and beyond the NREFT ROI, allowing for accurate characterization of ER leakage from β back-grounds into the NR signal region.
The NR response for LZ is only directly calibrated to approximately 80 keV nr1 using D-D with supplementary AmLi neutron calibrations.However, the NR models in NEST are based on a collection of all LXe light and charge yield measurements available in the existing literature, and the highest energy measurements included in NEST are the 330 keV yields from AmBe reported by Sorensen et al. [32].Therefore this analysis relies on extrapolating the NEST NR response beyond the in-situ LZ calibrations, specifically using recent ex-situ measurements of the NR response from D-T neutrons which provide information on the NR yields up to 426 keV [33].Ref. [33] suggests that the NEST models may be overestimating both the light and charge yields at higher energies beyond the D-D endpoint.Using those D-T measurements, the uncertainty in the NEST NR models beyond the in-situ calibration yields is constrained.However, when accounting for the reduction in total quanta, the change in the mean NR {log 10 (S1c), log 10 (S2c)} response is calculated to be <1%.The NEST models for the β ER and NR response and NR response uncertainty used in this analysis are compared to LZ calibration data and shown in Fig. 2. The impact from this uncertainty on WIMP sensitivity is tested by altering the NR response between -1σ and +3σ uncertainty.Because the discrimination between ERs and NRs significantly improves as energy increases, the NR response uncertainty has a negligible impact on the results of this analysis (<10% on the final WIMP sensitivity, for any combination of mass and operator).A series of data selection criteria are used to remove accidental coincidences of isolated S1s and S2s ("accidentals") from true single scatters at an efficiency greater than 99.5%.Sources of isolated S1 pulses are particle interactions in charge-insensitive regions of the TPC, Cherenkov and fluorescent light in detector materials, and dark-current pile-up.Isolated S2 pulses can be generated from radioactivity and electron emission from the cathode or gate electrodes, particle interactions in the gas phase or the liquid above the gate electrode, and delayed drift electron signals [34].The criteria for removing accidentals are tuned using side band events and are only applied to the search dataset after being finalized.
The accidental removal criteria are the same as the methods reported in Ref. [1], using relationships between pulse and event based quantities (such as drift time, the ratio between light collected in either the top or bottom PMT arrays, pulse width, the timing of PMT hits within the pulse, and hit pattern of the photons in the PMT arrays), and targeted individual sources of isolated S1 and S2 pulses by comparing to the expected single scatter behavior.The efficiency of the data selection criteria beyond the lower-energy ROI reported in Ref. [1] is evaluated using tritium and 220 Rn data for cuts targeting S1 pulses values and the combination of tritium and AmLi data for those targeting S2 pulses.
Following a similar approach as the LUX NREFT analysis [9], this analysis implements a boosted decision tree (BDT) to identify and remove γ-X events, which are the interactions of multiply scattering γ-rays classified as single scatters due to one or more scatters occurring in a region of the TPC from where charge cannot be collected [9].This leads to an attenuated S2 signal relative to the observed S1, increasing the ER leakage into the signal region.Two events are identified as γ-X and removed from the LZ NREFT search data, consistent with the model expectation of 1.6 events, all from the bottom of the FV.Fig. 3 shows the spatial distribution of the events passing all analysis criteria, highlighting the events removed by the γ-X classifier and the LXe Skin and OD veto systems.The NR acceptance of the BDT was calculated using simulated data to be 99.950±0.002%within the FV, and was validated using D-D and AmLi neutron calibration data.A fraction of events is classified as γ-X outside the FV, where the BDT performance is degraded due to the noisiness of S1 signals near the reflective TPC wall.Further details of the γ-X event topology, modeling procedure, and BDT results are discussed in Section III B 1.
The final efficiency of the data selection criteria, evaluated with 3 H, 220 Rn, and AmLi calibration data and including trigger and event reconstruction efficiency, is shown in Fig. 4 as a function of true NR energy.

III. MODELING
To simulate the background and signal components in the observable space, the BACCARAT package based on GEANT4 [35,36] is utilized, along with a bespoke simulation of the LZ detector response, which is fine-tuned using the NEST detector model.As part of this methodology, the uncertainties associated with the background components are included as constraint terms in a combined fit of the background model to the data.

A. Signals
As stated in Section I, the recoil spectra for each operator-mass combination are generated using a modified version of WimPyDD.This analysis considers both isoscalar and isovector interactions for WIMP masses between 10 and 4000 GeV/c 2 .As several operators have a dependence on the spin of the target nucleus, the response from each Xe isotope, weighted by its natural abundance, is incorporated into the overall response of the target.In this analysis, the isospin representation follows the methodology outlined by Anand et al. [14], previously used in Refs.[8,10] where c n and c p denote the coupling coefficient for interactions with neutrons and protons, respectively, and c s and c v represent isoscalar or isovector interactions, respectively.The normalization chosen here differs from some previous 1000 GeV, O s 13 FIG. 5. Example {log10(S1c), log10(S2c)} distributions for a 1000 GeV/c 2 WIMP-nucleon isoscalar interaction for the momentum-independent operator O1 (top), and momentumdependent operators O6 (middle) and O13 (bottom).Red and blue lines denote the flat ER and NR response regions as described in Fig. 2. Each pane shows the distribution for 100,000 WIMP nuclear recoils.
Each operator is considered as an independent interaction channel, where only contributions from that operator with no interference have been considered.In practice, this requires setting the a priori of the coefficient of the given operator to 1/m v , and all other coefficients to 0. The quantity m v is introduced to ensure the coefficients have dimensions of energy −2 and is set to the Higgs vacuum expectation value of 246.2 GeV/c 2 [7].This is a natural value for the normalization as it allows the coefficients to be expressed in terms of the Standard Model weak interaction mass scale.The normalization is re-quired due to the decision by Anand et al. to normalize the spinors by 4m N m χ , which allows for a dimensionless representation of the operators.By considering the case of a single operator O i dominating the interaction, the differential recoil rate scales with the square of c N i , where N ∈ {s, v}, such that where v min = q(m χ + m N )/(2m χ m N ) is the minimum velocity that produces a recoil energy E R , and is a form factor that depends on the nuclear physics.
A sample of the signal distributions in the {log 10 (S1c), log 10 (S2c)} observable space is shown in Fig. 5.

B. Backgrounds
The background model used in this analysis consists of 11 components.Tab.I lists the expected and fitted number of events for each component.Ref. [43] contains a complete discussion of the backgrounds in LZ for the data-taking period used for this analysis.Ref. [1] describes most of these background sources in detail, and the expectation for some backgrounds is mostly unchanged (such as 37 Ar, 8 B coherent neutrino-nuclear scattering, detector neutrons, and accidental coincidence of isolated S1s and S2s).The expected contributions from continuous ER sources -namely the "Flat ER" sum of 222 Rn, 220 Rn, and 85 Kr, as well as the ERs from solar neutrinos and double-β decays from 136 Xe -is increased from Ref. [1] due to the extended energy window of this analysis.Except for increasing the energy range, these models are unchanged from the previous analysis.The 127 Xe and 124 Xe models are expanded in this analysis, as the 127 Xe K-shell electron captures and 124 Xe KL, KM, and KN double electron capture signals are reconstructed partially within the search ROI 2 .
Two other sources of background are treated uniquely for this analysis: 125 I electron captures and Compton scattering γ-rays from trace levels of 40 K, 60 Co, 232 Th, and 238 U in the detector components [45] as well as 40 K, 232 Th, and 238 U from the cavern walls [46].
125 I is introduced into the TPC via neutron activation of 124 Xe into 125 Xe and its subsequent decay.With a 59.4 d half-life, 125 Xe produces unresolved multiple scatters with a 35.5 keV γ-ray in addition to an electron capture X-ray.The combined γ+L (40.4 keV) and γ+M (36.5 keV) decays contribute events into the search window of this analysis.The γ+K decay is outside of this search window and was used to infer the rate of 125 I.
Compton scatters from detector components were treated similarly to the Flat ER contributions in Ref. [1].However, the rate of these backgrounds increases with energy as higher energy decays have longer mean free paths in LXe [47].The longer mean free paths increase the probability of multiple scatters in the fiducial volume.Due to finite resolution, some of the events are not separable from single scatters, causing an increased rate of events to deviate from a standard ER single scatter response at higher energies.This is further obscured by the unique detector pathologies near detector components, such as poor light collection efficiency and field fringing, increasing the deviation.Therefore, detectorbased ERs are separated from other ER sources for the background fitting procedure in this analysis.These effects are typically subdominant in terms of ER leakage into the signal region since the ER and NR bands diverge at these energies unless a significant portion of the energy is deposited in a region from where the ionization signal cannot be collected.These γ-X events predominantly occur near the cathode and TPC walls.Because they are a unique background to high-energy searches and can be reconstructed near and below the NR signal region, they are considered separately in this analysis and described in the following section.
Figure 6 shows the {log 10 (S1c), log 10 (S2c} distribution of the 835 events which pass all selections, along with contours representing a 1000 GeV/c 2 O 6 isoscalar signal model (representative of signal models that peak at nonzero energy), and the background model.

TABLE I.
Expected and fitted numbers of events from the listed sources in the 60 d × 5.5 t exposure.The middle column shows the predicted number of events with uncertainties as described in the text.These uncertainties are used as constraints in a combined fit of the background model.The fit result is shown in the right column."Flat ER" represents the combination of 214 Pb, 212 Pb, and 85 Kr mixed in the LXe, while "Detector ER" represents electron recoils originating from radiogenic decays in detector materials.Both 37  FIG.6.The final high energy WIMP-search data after all cuts in {log10(S1c), log10(S2c)} space.The contours that enclose 1σ (dark) and 2.5σ (light) regions represent the following models: the shaded red region indicates the detector neutrons, the shaded orange region indicate the detector ERs, the blue region is the combined representation of all other ER models ( 214 Pb, 212 Pb, 85 Kr, 37 Ar, 125 I, 124 Xe, 127 Xe, 136 Xe, and ν ER) and the black dashed lines show a 1000 GeV/c 2 O6 isoscalar signal model.The solid red line corresponds to the NR median, while the red dotted lines represent the 10 -90% percentiles of the expected response.The model contours are produced with a linear scale for S1c prior to being cast into log and take into account all the efficiencies used in the analysis.Contours of constant recoil energy have been included as thin gray lines.Grayed regions at the left and top of the plot indicate parameter space outside the energy ROI.

γ-X Events
A γ-X event occurs when a γ-ray scatters multiple times in the TPC, but at least one scatter occurs in a region where electrons cannot be extracted (such as below the cathode electrode or in regions with significant electric field distortions).This leads to missing S2 pulses observed in the interaction, allowing these multiple scatters to be erroneously classified as single scatters with lower S2:S1 ratios than typical ER events, potentially mimicking an NR event.These pathologies can be ignored in a traditional SI WIMP search; the S2-prohibitive regions are typically spatially distant from the FV boundaries, so it becomes unlikely that a γ-X event can occur in the FV with an S1 signal small enough to be observed in a standard WIMP low energy ROI [1].However, in the extended energy window of an NREFT search, γ-X events become a significant NR-like background.
In LZ, the locations primarily responsible for the production of γ-X events are the reverse field region (RFR) below the cathode, and in close proximity to the walls of the TPC.Both regions have electric fields that do not allow the complete collection of charge.In the case of the RFR, electrons drift down instead of up.In the vicinity of the TPC wall, non-uniformities in the electric field direct drifting electrons to the wall, leading to the deple-tion of the charge signal.The most significant source of γ-X events below S1c = 600 phd is the electron capture decay of 127 Xe near the cathode, where the two-step deexcitation produces an X-ray photon (0.2-32.2 keV) and a γ-ray photon with energy O(100 keV).This same decay near the TPC wall and decays of rare-earth impurities in the cathode grid wires have minor contributions to the overall rate of γ-X events.
To generate γ-X events from 127 Xe and cathodic origins, we use a custom Monte Carlo simulation for γ-ray propagation through a realistic LZ geometry.By incorporating the measured activities of 127 Xe and radiogenic impurities in detector components, this simulation calculates the three-dimensional propagation of γ-rays from a given starting position, using the expected energydependent mean free path and Compton scattering crosssections in LXe [47,48].Light and charge yields are calculated using the photo-absorption and Compton scattering ER models from NEST.S1 position dependence and individual PMT channel areas are calculated using light collection maps generated with BACCARAT [35].Additionally, field magnitudes and electron drift paths throughout the TPC are calculated using maps generated from dedicated finite-element simulations [49], providing information about the regions along the cathode and wall where an S2 signal would escape detection.The simulations include the potential wells at the wall created by TPC field shaping rings that trap electrons, reducing the size of the S2 pulse.However, the rate of γ-X events at the wall is far lower than that at the cathode, given the much larger volume of the RFR.
A data-driven number of expected γ-X events was calculated using the γ-X model and the measured rate of 127 Xe, the dominant source of the background.The model was used to obtain the ratio of γ-X to the clean 203 keV peak from 127 Xe, which was then scaled by the size of the 203 keV peak in the observed dataset to yield an expectation of 1.6 ± 0.3 γ-X events in the ROI and FV.The rate of observed γ-X candidates rises in accordance with the model outside the FV as the edge of the TPC boundary is approached.
The γ-X simulations are used to train a multi-class BDT in order to classify events as ER, NR, 127 Xe RFR γ-X, 127 Xe wall γ-X, or cathode RFR γ-X.The BDT is provided with a seven-dimensional simulated dataset with the following features: -position-corrected S1 pulse area, -position-corrected S2 pulse area, -radial reconstructed position, -vertical position (drift time), -"cluster size", the dispersion of S1 light collected in the bottom PMT array, -"max channel fraction", the ratio of light observed in the brightest bottom array PMT to the total bottom array S1 area, -the ratio of light observed in the top and bottom PMT arrays.The final three S1 pulse-based features may be exploited to differentiate γ-X events because each of the multiple scatters contribute to the S1 pulse.None of the seven features individually can provide sufficient γ-X rejection power.The information contained in the correlations of the hit pattern features with pulse areas and event positions, however, may be harnessed by a multivariate tool such as a BDT for the clean removal of γ-X events.The accuracy of true single scatter and γ-X identification is assessed by cross-validating BDT predictions for ten non-overlapping datasets.In this procedure, the entire simulated dataset is split into ten equally-sized portions (80,000 events), and each portion is taken as the validation set for a BDT that is trained on the remaining nine portions.The parameters and classification thresholds of the ten BDTs are identical to those of the final BDT deployed on the EFT search data.The general performance on the validation datasets is summarized in Tab.III B 1 with errors given by the standard deviations of predicted counts across the ten BDTs.A high averaged single scatter acceptance is seen within the FV and ROI, which remains high even at the boundaries of the FV.To validate the high acceptance of simulated single scatters in observed data, the BDT is deployed on 220 Rn calibration data in the S1c < 1000 phd range.Under the most conservative assumption that all 25,000 events originating from the 220 Rn decay chain are true single scatters, the single scatter acceptance rate of the BDT is 99.92%.The BDT is also deployed on DD and AmLi neutron calibration data to verify NR acceptance.No DD events are removed by the BDT, and only AmLi events below the NR band that are indicative of multiply scattering neutrons were removed (multiply scattering neutron events have the same S1 hit pattern characteristics of γ-X events).Finally, the γ-X rejection rate in the 1000 phd < S1c < 2000 phd side band of the WIMP search data is tested.This side band has a significant number of γ-X events from the 375 keV γ-ray of the 127 Xe electron capture decay.The BDT removes 73% of the 180 γ-X events on and below the NR band.The 27% γ-X misclassification rate is attributed to differences in the data and the detector response model (used to train the BDT), which is not tuned or validated in this S1 range.
The assessment using calibration data shows that the BDT has an acceptance larger than 99.9% for ER and NR single scatter events throughout the FV.A lower bound on the γ-X rejection rate of 73% is obtained from a side band above the ROI, with the true rejection rate in the ROI expected to approach the value from simulations (99.6%) at lower energies due to the better match between data and the detector response model.Two events are removed by the BDT in the search dataset, consistent with the expectation of 1.6 γ-X events found using the measured activities and the rate at which they produce γ-X in the custom Monte Carlo simulation.This consistency, in addition to the high γ-X rejection rate of the BDT classifier, removes the need to incorporate a γ-X model in the fits to the observed data.

IV. RESULTS
No significant evidence of an excess is found among the NREFT operators in the isoscalar and isovector bases for both elastic and inelastic DM in any of the models tested.Comparisons of the reconstructed energy distributions between the data and the background model using unbinned Kolmogorov-Smirnov (Anderson-Darling) tests yield p values of 0.392 (0.25), showing consistency between the data and the background-only scenario.Upper limits on the DM coupling strengths for each NREFT interaction are presented in Section IV A for elastic scatters and Section IV B for DM upscattering to a heavier state.
The upper limits are obtained by defining an extended unbinned profile likelihood statistic in log 10 (S2c)-S1c space, which is used to construct two-sided bounds at the 90% confidence level following Refs.[1,37].The resulting coupling strength limits are cast in the dimensionless form (c N i × m 2 v ) 2 , where m v = 246.2GeV/c 2 is the Higgs vacuum expectation value [14].

A. Elastic
The upper limits for the coupling strengths of DM scattering elastically via the operators O 1,3−15 are shown in Fig. 7 and Fig. 8 for isoscalar and isovector interactions, respectively.A power constraint is applied on all operators (in the 17-30 GeV/c 2 mass range) to restrict the upper limit falling 1σ below the median expectation due to background fluctuations.The constrained limit corresponded to a critical alternate hypothesis power of π crit = 0.16 [1,37,50].Observed upper limits are compatible with background-only expectations to within 1σ for the majority of the operator-mass combinations.A few operator-mass combinations, e.g.O 3,13,15 above 30 GeV/c 2 , are found to have limits weaker than the median expectation.These discrepancies are within 2σ and are generally caused by ER background leakage into the NR region occupied by these momentum-suppressed operators with highly peaking spectra (see Fig. 1).The deviation from the ER band of the two most egregious outlier events is consistent with unresolved multiple scatter from detector-based ER decays.
Consistency with the first LZ result (Ref.[1]) is established with O 1 , the SI operator that couples solely to the total number of nucleons.Unaffected by the ⃗ v ⊥ and ⃗ q degrees of freedom, O 1 and O 4 yield some of the most stringent constraints on the couplings of nearly all the operators.The most stringent limit is set by O 1 and the limit on O 4 is only exceeded by those of O 11 and O 12 .The nuclear response of O 11 is similar to that of the SI operator O 1 at higher energies, however momentum-dependence causes a suppressed rate at low energies.O 12 is an example of an operator for which the rate is enhanced by positive parity nucleon velocity contributions that are summed over the composite nucleus [14].While, the exclusion curves for most operators have minima at WIMP masses of 30-50 GeV/c 2 , interactions that are suppressed by two powers of ⃗ q such as O 6 and O 15 attained minima in their coupling strengths at WIMP masses of 200-300 GeV/c 2 ; these are the operators that benefit the most from an extended energy window.
Results of the XENON100 [10], LUX [9], and Pan-daX [8] isoscalar analyses are also shown in Fig. 7 for comparison.The LUX constraints are only available for a WIMP mass of 1 TeV/c 2 since they originate from an inelastic DM search from which the data for zero mass splitting were used.All results are normalized to the dimensionless form (c N i × m 2 v ) 2 , and the differing normalization conventions among experiments are presented in Appendix A. Several other results are omitted in the comparison, primarily due to their choice of presenting limits in the proton and neutron bases instead of the isoscalar and isovector bases.Previous searches for NREFT interactions also used one-body nuclear density matrices for Xe that have since been updated [51].This analysis uses the updated matrices with the effect of altered event rates for some operators, notably a decrease in the event rate for O 13 , leading to an upper limit weaker than previous results in some cases.

B. Inelastic
Limits on the coupling strengths c I i for inelastic interactions are presented in Fig. 9 (isoscalar) and Fig. 10 (isovector) for values of the mass splitting δ spanning up to 250 keV and WIMP masses from 400 GeV/c 2 to 4 TeV/c 2 .For lighter WIMPs, an increasingly larger fraction of the scattering energy is required to transition to the heavier state.Therefore, the inelastic rates for lighter WIMPs increasingly fall below the LZ energy threshold, and so WIMPs lighter than 400 GeV/c 2 are not considered.The same kinematic suppression leads to weaker limits for larger values of δ at all WIMP masses [20].The observed inelastic upper limits do not drop below 1σ of the background-only expectations; therefore, no power constraint is required.

V. CONCLUSION
This article presents the results of a search for a wide range of dark matter scenarios using the first 60 liveday run of the LZ detector.The study extends the energy window of Ref. [1] to include nuclear recoils of up to around 270 keV nr (defined as where the acceptance starts to fall below 90%), which requires modeling and removal of the γ-X background.A frequentist statistical analysis tested the data against 14 WIMP-nucleon operators generated by an NREFT, for both elastic and inelastic scattering, and for isoscalar and isovector couplings to the xenon nucleus.A total of 56 distinct interactions were tested, corresponding to a set of nuclear recoil spectra that span the entire energy window used in the analysis.No significant evidence of an excess is observed for any model.
The coupling strengths of all possible elastic (inelastic) interactions between nuclei and dark matter with mass 10-10 4 GeV/c 2 (400-4000 GeV/c 2 ) are significantly constrained.LZ provides the strongest measured upper limits for nearly all the models tested.In particular, models with momentum-suppressed recoil spectra were tightly constrained due to the extended energy window that provided a higher efficiency for the resulting signal events.This result paves the way for NREFT interactions to be constrained further by future searches that leverage target nuclei with different isospin properties to xenon [52].

FIG. 2 .
FIG.2.Calibration events in {log10(S1c), log10(S2c)} space from tritium (light blue), D-D neutrons (orange), and 212 Pb (green) in and above the ROI.The mean ER and NR responses from NEST for flat energy spectra are shown in dark blue and red, respectively; the dashed lines indicate the 10% -90% percentiles of the expected response.Additionally, the region highlighted in pink denotes the shift in the 10% -90% percentiles when considering the ±1σ uncertainty in the NEST mean NR response beyond the D-D endpoint used for this study.Gray contours are lines of constant reconstructed energy, labeled for both ER and NR interactions.The gray dashed horizontal line denotes the upper log10(S2c) bound used in this analysis, however, the calibration events in this region were used to constrain the ER response model.

FIG. 3 .
FIG.3.Data in reconstructed r 2 and z after all analysis cuts.Black (gray) points show the data inside (outside) the fiducial volume after all cuts and vetoes have been applied.Red crosses and blue circles show events vetoed by a prompt Skin and OD signal, respectively.Solid green diamonds indicate the events removed by the γ-X BDT cut after all other cuts have been applied.Hollow diamonds indicate events outside the FV classified as γ-X.The solid line shows the fiducial volume definition, and the dashed line shows the extent of the active TPC.

FIG. 4 .
FIG.4.Signal efficiency as a function of the NR energy from the trigger (blue), the three-fold coincidence and 3 phd threshold on S1c (orange), the single-scatter (SS) reconstruction and analysis cuts (green), and the search ROI in S1c and S2c (black).The low energy behavior is shown in the inset, where the dotted line at 5.5 keVnr indicates the nuclear recoil energy at which the efficiencies equal 50%.The uncertainty on the detection efficiency (gray region) was assessed with 3 H, 220 Rn, and AmLi calibration data.

15 FIG. 10 .
FIG.10.The 90% confidence limit (solid lines) of the dimensionless isovector WIMP-nucleon couplings for each of the fourteen NREFT interactions.The dotted lines show the medians of the sensitivity projection and the shaded bands correspond to the 1σ sensitivity band.The upper limit is evaluated for WIMP masses of 400 GeV/c 2 , 1000 GeV/c 2 , and 4000 GeV/c 2 , for values of the mass splitting δ of 0 keV (purple), 50 keV (blue), 100 keV (green), 150 keV (yellow), 200 keV (orange), and 250 keV (red).
Ar and the detector neutrons have non-gaussian constraints and are totaled separately.Values with a fit result of zero are set to have no lower uncertainty.

TABLE II .
Confusion matrix after FV and ROI cuts, averaged over ten BDTs (each evaluated on a different dataset), showing the correct identification rate (%) across the two classes on the diagonal and the misclassification rate on the off-diagonals.